Patrice Bertail

Mohammed Bouchouia

Ons Jelassi

Jessica Tressou

Mélanie Zetlaoui
email: melanie.zetlaoui@gmail.com

Scaling by subsampling for big data, with applications to statistical learning

Keywords: Scaling, big data, Subsampling, Convergence rate estimation, Confidence intervals in statistical learning, Out-of sample error, EMNIST digits VeReMi

Handling large datasets and calculating complex statistics on huge datasets require important computing resources. Using subsampling methods to calculate statistics of interest on small samples is often used in practice to reduce computational complexity, for instance using the divide and conquer strategy. In this article, we recall some results on subsampling distributions and derive a precise rate of convergence for these quantities and the corresponding quantiles. We also develop some standardization techniques based on subsampling unstandardized statistics in the framework of large datasets. It is argued that using several subsampling distributions with different subsampling sizes brings a lot of information on the behavior of statistical learning procedures: subsampling allows to estimate the rate of convergence of different algorithms, to estimate the variability of complex statistics, to estimate confidence intervals for out-of-sample errors and interpolate their values at larger scales. These results are illustrated on simulations, but also on two important datasets, frequently analyzed in the statistical learning community, EMNIST (recognition of digits) and VeReMi (analysis of Network Vehicular Reference Misbehavior).

Introduction

The capacity to collect data has increased faster than our ability to analyze big datasets with a huge number of individuals. The standard statistical tools or statistical learning algorithms, like machine-learning procedures, maximum likelihood estimations, or general methods based on contrast minimization, are time-consuming in terms of optimization despite their polynomial complexities or require to access the data too many times. For these reasons, these approaches may be unsuitable for big data problems. To remedy the apparent intractability problem of learning from databases of explosive sizes and break the current computational barriers, we propose in this paper to use some variations of subsampling techniques studied in [START_REF] Politis | Subsampling[END_REF] and [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF]. Such approaches have been implemented in many applied problems and developed for instance in [START_REF] Kleiner | A scalable bootstrap for massive data[END_REF]. In the framework of big data, they are also at the core of some recent developments in survey sampling methods [START_REF] Clémençon | Scaling up m-estimation via sampling designs: The horvitz-thompson stochastic gradient descent[END_REF][START_REF] Bertail | Tail index estimation based on survey data[END_REF][START_REF] Bertail | Empirical Processes in Survey Sampling with (Conditional) Poisson Designs[END_REF] applied to statistical learning procedures.

The universal validity of the subsampling methods is proved in [START_REF] Politis | Large sample confidence regions based on subsamples under minimal assumptions[END_REF] and further developed in [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF][START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF] for general converging or diverging statistics. A good survey on the power of subsampling methodology is given in [START_REF] Politis | Subsampling[END_REF]. More precisely, the subsampling distribution constructed with a much smaller size than the original one is a correct approximation of the distribution of the statistic of interest (possibly with an unknown rate of convergence), if the latter has a non-degenerate distribution, that is continuous at some point of interest. As mentioned in [START_REF] Hall | A short prehistory of the bootstrap[END_REF], such methods, close to bootstrap and subsampling ideas, were already proposed in the works of Mahalanobis in the 1940s (see the reprint [START_REF] Mahalanobis | Recent experiments in statistical sampling in the indian statistical institute[END_REF] but were abandoned because of the cost of paper: at that time, calculations were carried out by hand. They have also been developed by Bretagnolle (1983) and [START_REF] Bickel | Richardson extrapolation and the bootstrap[END_REF]. See also the discussions about interpolations and extrapolations in [START_REF] Bertail | Extrapolation of subsampling distribution estimators: The i.i.d. and strong mixing cases[END_REF], [START_REF] Bertail | Second-order properties of an extrapolated bootstrap without replacement under weak assumptions[END_REF], when the computer capacities were not sufficient to handle even moderate sample sizes. Such methods are themselves related to well-known numerical methods (see for instance [START_REF] Isaacson | Analysis of numerical methods[END_REF][START_REF] Har-Peled | Geometric Approximation Algorithms[END_REF].

Most of these methods based on subsampling rely on an adequate standardization of the statistics of interest. Such standardization may be hard to obtain for complicated procedures including statistical learning procedures. It is even more complicated to extrapolate to very large sample sizes. Indeed, the extrapolation of the distribution of a statistic from smaller scales to a large one requires the knowledge of the rate of convergence τ n of the procedure of interest or at least an estimator of the latter, with n being the size of the dataset. In many situations, this task is difficult because the rate itself depends on the true generating process of the data.

In this paper, we present a variant of the subsampling distribution estimation methodology studied by [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF][START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF] to derive a consistent estimator of the rate τ . and study its rate of convergence. We prove the asymptotic validity of the method to construct asymptotically valid confidence intervals and give some precise rate of convergence, assuming that the centering of the subsampling distribution satisfies some concentration inequality. What differs from previous works is the way the subsampling distribution of the statistic of interest is centered. It allows precise control of the rate of the approximation. In [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF][START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF], the centering quantity was the statistic of interest computed on the whole data set (which may not be achievable in practice for big data, when the size n is very large), and the associated rate of convergence was of the order of 1/ log n. Here, the centering quantity is computed as the mean of the statistics of interest obtained on subsamples of size b n . This yields for subsampling distribution, the expected rate of convergence of order b n /n plus a bias term of order (in regular cases) 1/ √ b n which is the rate at which the true distribution based on a sample of size b n approaches the asymptotic distribution. In that case, we obtain an optimal rate of n -1/4 for the choice b n = n -1/2 . Moreover, we show that the mean of the subsampling distribution (eventually approximated by Monte Carlo simulations) yields attractive hyper-efficiency properties for some slow algorithms (with a rate slower than √ n). This estimator satisfies precisely the required concentration inequality. We then use the extrapolation of this estimation for large datasets to construct confidence intervals for many procedures that are difficult to analyze otherwise. We show how this idea may be used in the case of ma-chine learning procedures to obtain confidence intervals for general risks. By enabling the construction of confidence intervals for the test error rates or out-of-sample error rates of various algorithms, it becomes easier to compare them. We also show how it is possible to practically integrate the dynamic aspect of the big data environments (especially in the case of streaming data flows).

These subsampling techniques are then implemented with potentially timeconsuming procedures (k-nearest neighbors, random forest, neural nets,...) first on simulated data and next on two real databases, the EMNIST dataset, and the VeReMi dataset. The implementation of simulated data is performed with the software R on a standard machine and the implementation of real data illustrations with Python on an Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz machine with 144 cores and 250GB RAM. We estimate the rate of convergence of several algorithms and obtain confidence intervals for the out-of-sample risks of several standard algorithms. An interesting byproduct of the study is that using different subsampling sizes also allows for detecting the instability of the procedures.

The article is organized as follows. Section 2 presents the state of the art of the subsampling methods and introduces some estimators of the convergence rate of the sample statistic distribution. Section 3 presents the two main results and a discussion on the choice of the subsampling sizes. First, we prove the asymptotic validity of the method to construct asymptotically valid confidence intervals. Then, we also prove that the mean of the randomized subsampling distribution yields attractive hyper-efficiency properties for some slow algorithms. Section 4 presents applications on simulated data and on the VeReMi and EMNIST-digits datasets. Finally, the proofs are deferred to the Appendix section.

Subsampling methods for big data

Definition

In [START_REF] Politis | Large sample confidence regions based on subsamples under minimal assumptions[END_REF] (see also [START_REF] Politis | Subsampling[END_REF] for more a developed framework and references) a general subsampling methodology has been proposed for the construction of large-sample confidence regions for an unknown parameter θ = θ(P ) ∈ R q under very minimal conditions. Considering X n = (X 1 , ..., X n ) an i.i.d. sample, the construction of confidence intervals for θ requires an approximation to the sampling distribution under P , generally unknown, of a standardized statistic T n = T n (X n ). This statistic is assumed to be consistent for θ at some known rate τ n . For example, in the statistical learning methodology, θ may be the Bayes Risk and T n the estimated risk linked to a given algorithm (see section 3.3 for illustration). In the framework of prediction, θ may be a value to predict and T n a predictor.

To fix some notations, assume that there is a non-degenerate asymptotic distribution for the centered re-normalized statistic τ n (T n -θ), denoted by K(x, P ), continuous in x, such that for any real number x,

K n (x, P ) ≡ Pr P {τ n (T n -θ) ≤ x} -→ n→∞ K(x, P ). (A1)
Recall that, in the framework of resampling techniques, we are interested in estimating the distribution of the original statistics K n (x, P ) or more generally the sequence of distributions K m (x, P ), m ≤ n. Indeed for big data, the statistic T n itself may be difficult to compute on the whole database (see the applications to the EMNIST and reduced EMNIST dataset in part 4.2). Then the subsampling distribution with subsampling size b n , is defined by

K bn (x | X n , τ . ) ≡ q -1 q i=1 1{τ bn (T bn,i -T n ) ≤ x}, (1)
where q = n bn and T bn,i is a value of the statistic of interest, calculated on a subset of size b n chosen from X n . It was shown in [START_REF] Politis | Large sample confidence regions based on subsamples under minimal assumptions[END_REF] that the the subsampling methodology is asymptotically valid. Precisely, under the following assumptions on b

n b n -→ n→∞ ∞, b n n -→ n→∞ 0, ( A2
)
and

τ bn τ n -→ n→∞ 0. (A3)
we have

K bn (x | X n , τ . ) -K n (x, P ) -→ n→∞ 0,
uniformly in x over neighborhoods of continuity points of K(x, P ). The same result holds if n is replaced by any sequence m n and the corresponding subsampling size b n satisfies conditions (A2) and (A3).

Remark(On conditions A2-A3)

The key point for this result is based on the fact that, when T n is replaced by θ in equation ( 1), one obtains a U -statistic of degree b n whose variance is of order bn n . Then the condition A2 ensures that the mean of this U -statistic K bn (x, P ) converges to a limiting distribution and that the variance of the U -statistic converges to 0. The condition A3 allows replacing the true value of the parameter by T n . In that case, A3 ensures that this replacement does not affect the limiting distribution. When choosing an adequate re-centering that may differ from T n (for instance the mean or median of the subsampling distribution), the condition A3 may be completely dropped as discussed below.

For large databases, computing q = n bn values of the statistics T bn,i may be unfeasible. In this case, it is recommended to use its Monte-Carlo approximation

K (B) bn (x | X n , τ . ) = B -1 B j=1 1{τ bn (T bn,j -T n ) ≤ x}, (2)
where {T bn,j } j=1,...,B are now the values of the statistic calculated on B subsamples of size b n taken without replacement from the original sample. It can be easily shown by incomplete U -statistics arguments that if B is large then the error induced by the Monte-Carlo step on the subsampling distribution is only of size

K (B) bn (x | X n , τ . ) -K bn (x | X n , τ . ) = O P 1 √ B ,
where the notation O P (.) refers to stochastic boundedness. We recall that Z B = O P (a B ) means that, for any ε > 0, there exists a finite δ ε > 0 and a finite

N ε > 0 such that, ∀B > N ε , Pr(|Z B /a B | > δ ε ) < ε.
Then, if the error of K bn (x | X n , τ . ) on the true distribution is controlled, it is always possible to find a value of B (eventually linked to n) such that the Monte-Carlo approximation is negligible. MacDiarmid's inequality applied to the indicator functions yields a precise control of the error at any probability error level δ.

The subsampling method is generally based on the centering by T n computed on the whole database. This centering may not be adapted for big data since the calculation of T n itself may be too complicated: the exact size may be unknown or the complexity of the algorithm and the cost induced by retrieving all the information may be too high. In the subsampling method, the main reason for using the centering by T n is simply due to the fact that under the condition A3, the convergence rate of T n , τ n , is faster than τ bn . Indeed :

τ bn (T bn,j -T n ) = τ bn (T bn,j -θ) + τ bn (θ -T n ) = τ bn (T bn,j -θ) + O P τ bn τ n = τ bn (T bn,j -θ) + o P (1).
This suggests using any centering whose convergence rate is faster than τ bn . This is, for instance, the case if one constructs a subsampling distribution without any centering or standardization with a subsampling size m n >> b n such that bn mn → 0 and τb n τm n → 0. In this case we have that 1 B B j=1 T mn,j which is a proxy of 1 q q j=1 T mn,j (with an error of size 1/ √ B) converges to θ at a rate at least as fast as τ mn (provided that the expectation of these quantities exists). The same results hold if one chooses the median rather than the mean of the T mn,j s (when considering the mean this amounts to recenter at the median of means) as considered in [START_REF] Laforgue | On medians of (Randomized) pairwise means[END_REF].

Rate of convergence for subsampling distributions

To our knowledge, the precise rate of convergence of subsampling distributions and their Monte-Carlo approximations has not been investigated except in some very specific cases (and especially the mean). Actually, in general, this requires more precise control of the centering factor and of the modulus of continuity of the asymptotic distribution as well as some control on the rate of approximation to the asymptotic distribution.

In the following, we denote by θ n any centering such that

τ bn (θ -θ n ) = o P (1). ( 3
)
For reasons that appear clearly in the proofs and for further applications in statistical learning, we assume that there is a concentration inequality for this estimator of the following form: for some η > 0, there exists some universal constants

C i > 0 , i = 1, 2 such that Pr τ bn | θ n -θ| > η ≤ C 1 exp(- n b n C 2 η 2 ). (A4)
In general, because the right rate of convergence of θ n is much more rapid than τ bn , such concentration inequality holds for regular statistics. This is the case of the mean or the moments for instance, under the existence of some exponential moments. We will show later that the mean of the statistics computed on (all or a sufficiently large number of) subsamples of well-chosen size b n satisfies such a requirement under minimal variance assumptions (see [START_REF] Laforgue | On medians of (Randomized) pairwise means[END_REF] for other concentration results for the median of means).

As in [START_REF] Götze | Adaptive choice of bootstrap sample sizes[END_REF], we need to control the (deterministic) convergence rate of the true distribution to the asymptotic distribution (this actually plays the role of the bias in approximating the true distribution). In general, this rate is given by some Berry-Esseen theorems or Edgeworth expansions. For this we assume that, for any n large enough, uniformly in x,

K n (x, P ) -K(x, P ) = O 1 n β , for some β > 0. (A5)
In regular cases, Berry-Esseen theorems yield typically a rate of approximation n -1/2 with β = 1/2 but for instance, for symmetric statistics (or asymptotically χ 2 distribution), we rather expect β = 1 .

We finally assume that we can locally control the modulus of continuity of the asymptotic distribution at points of continuity K(., P ), by some increasing function null at 0. To simplify we assume that the distribution is Lipschitz in neighborhoods of continuity points (similar results but with different rates can be obtained under Hölder assumptions). There exists L > 0, such that, at any point x of continuity of K(., P ) and any y in a neighborhood of x.

|K(y,

P ) -K(x, P )| ≤ L|y -x|. (A6)
Since most of the time, the asymptotic distribution is differentiable with a bounded derivative, the distribution is Lipschitz and this hypothesis is trivially satisfied.

For simplicity, we use the same notation as before and now define the subsampling distribution with the centering θ n as

K bn (x | X n , τ . ) ≡ q -1 q i=1 1{τ bn (T bn,i -θ n ) ≤ x},
and its Monte-Carlo approximation

K (B) bn (x | X n , τ . ) = B -1 B j=1 1{τ bn (T bn,j -θ n ) ≤ x}.
The following result gives a new rate of convergence for K bn (x | X n , τ . ) and K (B) bn (x | X n , τ . ) for the centering θ n . This rate will be denoted by

δ β (n) = b n n + 1 b β n .
Theorem 1 Assume that conditions A1 to A6 hold, then we have uniformly over the neighborhood of points of continuity of K (x , P ) (uniformly over R if K (x , P ) is continuous)

K bn (x | X n , τ . ) -K bn (x, P ) = O P b n n , K (B) bn (x | X n , τ . ) -K bn (x, P ) = O P b n n + O P 1 √ B .
Moreover, we have

K bn (x | X n , τ . ) -K(x, P ) = O P (δ β (n)) ,
and

K bn (x | X n , τ . ) -K n (x, P ) = O P (δ β (n)) ;
if in addition n bn = 0(B), then the above hold upon replacing K bn (x | X n , τ . ) with its Monte-Carlo approximation

K (B) bn (x | X n , τ .
). This result gives some precise results on how to choose the subsampling size and the number of replications by optimizing the rate of the approximation. Typically for β = 1/2, we can choose the optimal value b n = n 1/2 and B = n 1/2 which actually will drastically reduce the computation costs in comparison to bootstrap procedures. In this case, δ β (n) = n -1/4 . For β = 1 (the symmetric case), we can choose b n = n 1/3 and B = n 2/3 , and in this case, we get a better rate δ β (n) = n -1/3 . It can be seen that the smaller the bias in A5, the better the approximation, but in this case, the Monte Carlo step will require more computations.

Estimating the convergence rate

The main drawback of this approach is the knowledge of the standardization (or rate) τ n . However, this rate may be easily estimated at least when the rate of convergence is of the form τ n = n α L(n) as shown in [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF]). Here α is an unknown real and L is a normalized slowly varying function, that is, such that L(1) = 1 and for any λ > 0, lim x→∞ L(λx) L(x) = 1 (see [START_REF] Bingham | Regular variation[END_REF]). For simplicity, we will now assume that τ n = n α . The general case τ n = n α L(n) may be treated similarly with additional assumptions on the slowly varying function (see [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF]). In any case, the estimator proposed in [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF] may be used in our framework. We now propose a simplified approach that is attractive in practice.

First, construct the subsampling without any standardization and denote for simplicity

K bn (x | X n ) ≡ K bn (x | X n , 1)
the subsampling distribution of the root (T n -θ). Then we have

K bn (x τ -1 bn | X n ) = K bn (x | X n , τ . ) (4)
Let denote now F -1 (t) the quantile transformation, i.e., F -1 (t) = inf {x : F (x) ≥ t} for a given distribution F on the real line and a number t ∈ (0, 1). The following Lemma extends Lemma 1 of [START_REF] Bertail | Extrapolation of subsampling distribution estimators: The i.i.d. and strong mixing cases[END_REF] by providing a rate of convergence for quantiles of subsampling distribution under natural assumptions on the limiting distribution.

Lemma 1. Assume that K(x, P ) is strictly increasing at least in the neighborhood of the quantile of interest K -1 (t, P ) . Then under the conditions of Theorem 1, we have

K -1 bn (t | X n , τ . ) = K -1 (t, P ) + O P (δ β (n) -1 ).
Now, it is easy to see with the rate of convergence obtained in Theorem 1 that we have by Lemma 1,

K -1 bn (t | X n , τ . ) = τ bn K -1 bn (t | X n ) (5) = K -1 (t, P ) (1 + O P (δ β (n))) , (6)
yielding log |K -1 bn (t | X n )| = log |K -1 (t, P )| -α log(b n ) + O P (δ β (n)) .
If two different subsampling sizes satisfy the conditions stated before and are such that b

n1 = b n , b n1 /b n2 = e, then one gets log |K -1 bn 1 (t | X n )| -log |K -1 bn 2 (t | X n )| = α + O P (δ β (n)) ,
uniformly in a neighborhood of t. Using sample sizes of the same order avoids the complicated constructions used in [START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF] and suggests that the parameter α may be simply estimated by averaging this quantity over several subsampling distributions. The rate that we obtain here is clearly better than the one obtained in [START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF] : this is due to the concentration hypothesis (A4) and hypothesis (A6). Again for the regular case β = 1/2, we can choose the optimal value b n1 = n 1/2 and obtain a rate for α equal to n -1/4 which clearly improves over the log(n) -1 rate obtained in [START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF]. Due to the "bias" term (involved by the asymptotic approximation rate in A5), this rate can not be improved.

Computing these two subsampling distributions mainly requires the computation of B × b n1 (1 + e) values of the statistic of interest (whose computation may be easily distributed). Thus the resampling size should be chosen large enough not to perturb too much the subsampling distributions but sufficiently small so that the cost of computing these quantities is small in comparison to the global cost of computing a single statistic over the whole database. Similarly as before, if we choose B = n 1/2 then we may replace the true subsampling distribution with the Monte-Carlo approximation without changing the rate. However, for some very large datasets, one may wish to choose a smaller B to have a lower computation cost, at the price of a worse approximation.

Notice that the previously described estimator of the rate parameter α involves the difference of two log-subsample quantiles, which has the drawback of depending on the subsample centering. Such centering can be removed by considering an estimator that involves the differencing of subsample quantiles before applying logarithms. For any 0 < t 1 < 1/2 < t 2 < 1, at which the asymptotic distribution is strictly increasing, we have

log K -1 bn (t 2 | X n ) -K -1 bn (t 1 | X n ) = log K -1 (t 2 , P ) -K -1 (t 1 , P ) -α log(b n ) + O P (δ β (n)) . For b n1 = b n → ∞ and b n2 = b n1 /e, it follows that log K -1 bn 1 (t 2 | X n ) -K -1 bn 1 (t 1 | X n ) K -1 bn 2 (t 2 | X n ) -K -1 bn 2 (t 1 | X n ) = α + O P (δ β (n)) . (7)
By looking simply at two subsampling distributions, it is thus possible to estimate the parameter α at a rate which is at least O P (δ β (n)). One may choose for instance t 1 = 0.75 and t 2 = 0.25 , corresponding to the log of interquartiles.

3. Subsampling with estimated rates : some rate of convergence

Confidence intervals based on subsampling

For a given estimator of τ n , typically τ n = n α , we will use

K n (x, P ) = Pr P { τ n (T n -θ) ≤ x}.
Theorem 2 Assume A1 holds for τ n = n α , for some α > 0 and some K(x, P ) continuous in x; assume also that assumption A2 holds. Let α = α + o P ((log n) -1 ), and τ n = n α. Then

∆ n = sup x |K bn (x | X n , τ . ) -K n (x , P )| = o P (1 ).
Let γ ∈ (0, 1), and let

c n (1 -γ) = K -1 bn (1 -γ | X n , τ . ) be the (1 -γ) th quantile of the subsampling distribution K bn (x | X n , τ . ). Then Pr P {τ n (T n -θ) ≥ c n (1 -γ)} -→ n→∞ γ. ( 8
)
Thus with an asymptotic coverage probability of 1 -γ, we have

-τ -1 n c n (1 -γ) ≤ θ -T n
and by symmetry,

θ -T n ≤ -τ -1 n c n (γ).
If in addition conditions A1-A6 hold, if we choose an estimator of α which satisfies

α = α + O P δ β (n) -1
then the rate of the subsampling approximation with the estimated rate becomes

∆ n = O P log(n)δ β (n) -1 .
Remark: Notice that estimating the rate as we did before, only results in a loss of log(n) in the subsampling distribution with an estimated rate. For β = 1/2 corresponding to smooth parameters, the optimal rate will be log(n)/n 1/4 . For moderate sample sizes, the subsampling approximation in Theorem 2 may be unsatisfactory. But, for very large datasets, it can still lead to an acceptable approximation, as will be seen in our applications.

Recall that

K -1 bn (1 -γ | X n , τ .
) is the (1 -γ) th quantile of the rescaled subsampling distribution. Just like in [START_REF] Hall | On the Number of Bootstrap Simulations Required to Construct a Confidence Interval[END_REF] assume that B is such that (B + 1) × γ is an integer (thus if γ = 5% or γ = 1%, B = 999 or B = 9999 is fine), the (1 -γ) th quantile is defined uniquely and equal to τ bn (T

((B+1)(1-γ)) bn -θ n ) where T ((B+1)(1-γ)) bn
is the (B + 1)(1 -γ) largest value over the B subsampled values. It then follows that the lower bound for θ -T n is given by

θ -T n ≥ - τ bn τn (T ((B+1)(1-γ)) bn -T n ). (9)
A straightforward application of this result is to compare the generalization capability of statistical learning algorithms (see section 3.3) when n is so large that the estimated risk T n of most algorithms, even with polynomial complexity, may be hardly computed in a reasonable time.

However, in some cases, it may be difficult to compute the statistics T n on the whole database. The preceding result also allows to build confidence intervals for θ based on a single subsampling realization say θ mn = T mn (X i1 , ..., X im n ) based on a different size m n such that n >> m n >> b n . For this, assume that (B + 1) × γ/2 is an integer. In that case, by combining (8) and ( 9), a two-sided confidence interval for θ is simply given by

θ mn - τ bn τmn T ((B+1)(1-γ/2)) bn -θ mn ≤ θ ≤ θ mn - τ bn τmn T ((B+1)γ/2) bn -θ mn .
This interval suffers of course from lower precision due to information loss in not being based on the whole data set. However, sometimes we can find an estimator θ n that can be computed on the whole database. For instance, the mean of a given subsampling distribution may be an adequate candidate in some cases (see part 3.2) with a rate that satisfies conditions (3) and A4. In that case, a rescaled confidence interval is simply given by

θ n - τ bn τn θ ((B+1)(1-γ/2)) bn -θ n ≤ θ ≤ θ n - τ bn τn θ ((B+1)γ/2) bn -θ n .
In our simulation studies and our applications, regarding the estimation of risks of a learning algorithm ( θn will be itself an estimated risk), the variability of the data may be so high that somehow there is very little difference between confidence intervals computed with T n or a subsampling estimator θ mn . To explain this phenomenon, recall that in statistical learning, the generalization error (see below) of an algorithm (or of an estimation procedure) can be decomposed into three terms

• the square of the bias, linked to erroneous assumptions in the learning algorithm (for instance a linear assumption while the model is not linear); • the variance, measuring the variability of the procedure or underlying estimator;

• the irreducible error of the model, that is, the variance of the intrinsic error.

From a statistical point of view, just consider the target model Y i = f (X i ) + i for some class of functions f ∈ F. The total L 2 -loss (or generalization error) of a predictor f (x) at a point x is

E(Y -f (x)) 2 = (E f (x) -f (x)) 2 + V ar( f (x)) + V ar( )
For big data, this last error is generally big due to the great variability of the observations whereas, at the same time, the bias and variance of the estimation procedure are small due to the large dataset.

Improving rates by means of subsampling estimators

In the particular case when the recentering estimator θ n is chosen to be the mean of the subsampling distribution, we show in the following theorem that there is a concentration phenomenon at a speed that is sufficient to ensure conditions (3) and A4. This result may be seen as a generalization of the hyper-efficiency results for pooled estimators computed on partitioned data obtained in [START_REF] Banerjee | Divide and conquer in nonstandard problems and the super-efficiency phenomenon[END_REF]. For slow procedures (i.e., procedures with a convergence rate slower than √ n), we prove this super-efficiency phenomenon that is, a convergence rate faster than the original rate and close to √ n, for the mean of subsampling distribution under a simple variance condition.

For simplicity, define the mean estimator θ q n = q -1 T bn,i

1≤i≤q
where q is either n bn or some deterministic B. In that case, the B subsamples are chosen uniformly over all possible subsamples of size b n .

Notice that, when we consider all subsamples, θ q n is nothing else than a U-statistic with a kernel T bn (.) of degree b n . Recall that if T bn is bounded by some constant M (which will the case in the statistical learning procedures that we will study later) and if the variance exists say V ar(τ bn T bn ) ≤ C < ∞, then Hoeffding proved a Bernstein type inequality (see also [START_REF] Arcones | A Bernstein-type inequality for U-statistics and U-processes[END_REF]) which becomes here

Pr(| θ q n -E(T bn )| > ε) ≤ 2 exp - n bn ε 2 M 2 2C/τ 2 bn + 2 3 M ε yielding (by changing ε into ε bn nτ 2 bn ) Pr( n b n τ bn | θ q n -E(T bn )| > ε) ≤ 2 exp   - ε 2 M 2 2C + 2 3 M ετ bn bn n  
Thus if b n is chosen such that τ bn bn n is bounded (which will be the case for the choice of b n controlling the bias and is always true when b n is chosen very small), then we will get an Hoeffding bound of type (A4). Moreover by a straightforward inversion of the Bernstein inequality, we have

| θ q n -E(T bn )| = O P ( b n n 1 τ bn + b n n )
which can be better than the rate of convergence of the original statistics τ n . We state this result under an unbiased condition (to avoid lengthy discussions on the form of the bias).

Theorem 3 Assume that the statistics of interest T n is an unbiased bounded estimator of θ has rate τ n ≤ √ n, under the assumptions A1-A2 and assuming that the asymptotic variance is bounded uniformly in n, that is,

V ar(τ n T n ) ≤ C.
Then the rate of convergence of θ q n , for q = n bn or q = B with B ≥ n/b n , is at least τ bn n bn . Moreover, this estimator satisfies the concentration inequality A4.

Remark In particular, if τ n = n 1/2 then the subsampling mean estimator θ q n has the same rate n 1/2 . But if τ n = n α , α < 1/2 then this estimator has a better rate of convergence given by n 1/2 /b 1/2-α n . We can even choose b n = log(n) and thus a rate close to √ n asymptotically. Of course, there is no free lunch. In that case, we need more computations of the subsampling estimator (at least n/ log(n) which may be too much in practice). Moreover, it is emphasized in [START_REF] Banerjee | Divide and conquer in nonstandard problems and the super-efficiency phenomenon[END_REF] that this estimator may not be locally regular which may be a drawback for some applications where uniformity is needed.

Some subsampling results tailored to statistical learning

Subsampling prediction error

Let D n = {(X i , Y i ), i = 1, ..., n} be i.i.d. with distribution P defined on (Ω, A, P), taking their values in some measurable product space X ×Y. The (X i , Y i )'s correspond to independent copies of a generic r.v. (X, Y ). A predictor is a measurable function φ : X → Y, x → y = φ(x). To measure the risk of a predictor, we introduce the loss function L : Y 2 → R which is assumed to be bounded by some constant M . The Bayes classifier is the one obtained by minimizing the expectation of the loss function over all classifiers :

φ * = arg min φ∈F E P L(Y, φ(X)).
However, in the estimation procedure, we will only minimize over a given class of function F corresponding to a specific algorithm. Moreover, since P is unknown, we will approximate the expected loss by the empirical one and consider the estimator φ n defined by

φ n = def φ n (F) = arg min φ∈F 1 n n i=1 L(Y i , φ(X i )).
One goal of statistical learning is to evaluate the generalization capability of the algorithm measured by the discrepancy between the optimal risk θ * = E P L(Y, φ * (X)) and the one evaluated on the resulting predictor φ n , say

∆ n = E P (L(Y, φ n (X))|D n ) -E P L(Y, φ * (X)).
Constructing confidence intervals for this quantity called prediction error can allow us to distinguish between different algorithms. In the following, it is assumed that ∆ n (recentered thus at 0) converges asymptotically to a distribution K(x, P ) at a rate τ n = τ n (P ), which is clearly unknown in most situations (even when one is able to obtain concentration inequalities). Notice that ∆ n which plays the role of T n in our first part is not really a statistics in the usual sense since we have taken expectation over (X, Y ). However, using subsampling ideas we can still get a subsampling approximation of the distribution of this quantity Following the subsampling ideas exposed before, for any subsampling set of size b n , D

(j) bn = {(X i , Y i ), i ∈ s (j) bn }, with s (j)
bn ⊂ {1, ..., n} for j = 1, .., q, we define the subsampling counterpart of

E P L(Y, φ n (X)|D n ) by E (j) bn = E P (L(Y, φ (j) bn (X))|D (j) bn ), evaluated at the estimator φ (j) b n = arg min φ∈F 1 b n i∈s (j) bn L(Y i , φ(X i )). Since E (j)
bn depends on the true distribution, we will estimate it by its empirical version computed on the set D

(j) bn = {(X i , Y i ), i ∈ s (j) bn } E (j) bn = 1 n -b n i∈s (j) bn L(Y i , φ (j) bn (X i )). Let's now define θ q n = Mean 1≤j≤q ( E (j)
bn ) and the subsampling distribution of the risk, on all q = n bn subsamples,

K bn (x | X n , τ . ) = q -1 q j=1 1{τ bn ( E (j) bn -θ q n ) ≤ x}.
As before, we introduce the approximate subsampling distribution based only on B simulations where now θ

B n = Mean 1≤j≤B ( E (j)
bn ) defined by

K (B) bn (x | X n , τ . ) = B -1 B j=1 1{τ bn ( E (j) bn -θ B n ) ≤ x},
where the E (j)

bn , j = 1, ..., B are taken at random uniformly on the set of all subsamples. We expect K (B) bn (x | X n , τ . ) to be an estimator of

K bn (x) = Pr τ bn (E (j) bn -θ * ) ≤ x ,
which is itself asymptotically close to the distribution of Pr P (τ n ∆ n ≤ x).

Then, we apply the same rate estimation procedure as before to compute an estimator of the convergence rate τ . , say τ . . Applying the same arguments as in Theorem 2 yields the following result.

Corollary 1 Assume that:

K (B) bn (x | X n , τ . ) Pr -→ n→∞ K(x, P ).
Moreover, under the same hypotheses as in Theorem 2, we have, with an estimated rate of convergence,

K (B) bn (x | X n , τ . ) Pr -→ n→∞ K(x, P ) yielding a confidence interval for ∆ n of level 1 -γ given by τ -1 n c n (γ/2) ≤ ∆ n ≤ τ -1 n c n (1 -γ/2)
where c n (t) is the quantile of order t of the distribution

K (B) bn (x | X n , τ . ).
Remark By the same arguments as in the proof of corollary 1,

θ q n = Mean 1≤i≤q ( E (i)
bn ) satisfies a concentration inequality around the true risk, provided that the variance of the risk estimator has an asymptotic variance uniformly bounded. This is automatically satisfied since the loss function L is assumed to be bounded. The problem reduces to obtaining a concentration result for

K (B) bn (x ) = B -1 B j=1 1{τ bn (E (j) bn -θ * ) ≤ x},
which is simply a U -statistic of degree b n .

Subsampling out-of-sample error or generalization error

It is known in statistical theory that conditional risk or measurement of prediction error may be too optimistic and could lead to the choice of an algorithm that overfits the data. For this reason, it is often recommended to split the data (of size n) into two parts, the training set called D n,T r of size n T rain = p 0 n, for some p 0 > 1/2 and the testing set say D n,T est of size n T est = (1 -p 0 )n. To avoid unnecessary problems and truncation, we assume that these sizes are integers. In our simulations we will choose p 0 = 0.7. The training and test sets are generally selected at random (practically using first a random permutation of the data and selecting the n T rain first units for the training set). In this framework, we are interested in the out-sample error estimator or generalization error given by

T n = 1 (1 -p 0 )n j∈Dn,T est L(Y j , φ nT rain (X j )))
where now φ nT rain is determined on the training set D n,T r that is

φ nT rain = arg min φ∈F 1 p 0 n i∈Dn,T r L(Y i , φ(X i )).
T n may be interpreted as an estimator of the predictive performance of the algorithm. For a convergent algorithm, that is, an algorithm which yields a loss asymptotically close to the Bayesian risk, we expect this quantity to converge to the unconditional full risk θ = E P L(Y i , φ * (X i )). In this framework, assuming that T n -θ is converging asymptotically to a non-degenerate distribution, we can apply straightforwardly the results of parts 2 and 3. Indeed if we choose a subsample D

T (k) bn = 1 (1 -p 0 )b n j∈D (k) bn ,T est L(Y j , φ (k) bnT rain (X j )))
where now φ (k) bnT rain is given by

φ nT rain = arg min φ∈F 1 p 0 n i∈Db n,T r L(Y i , φ(X i )).
The subsampling distribution is then the usual one and we can obtain confidence intervals for θ but also improve the convergence rate of T n when it is not converging at rate √ n.

How to choose the "optimal" subsampling sizes

The choice of the subsampling size is a delicate subject that has been discussed in very few papers including [START_REF] Bickel | On the choice of m in the m out of n bootstrap and confidence bounds for extrema[END_REF]; [START_REF] Götze | Adaptive choice of bootstrap sample sizes[END_REF]; [START_REF] Bickel | Subsampling methods for genomic inference[END_REF]; [START_REF] Politis | Subsampling[END_REF] Chap.9. Our preceding results were essentially asymptotic. For instance, for β = 1/2 the optimal choice is of the form b n = C √ n for some constant C. But in practice the choice of the constant is crucial... The main idea underlying most propositions is to construct several subsampling distributions by using two different subsampling sizes say b n and kb n for k ∈]0, 1[ (we recommend due to our preceding results k = 1/e). It is easy to see that when the subsampling distribution is a convergent estimator of the true distribution then the distance d between the subsampling distribution and the true one is stochastically equivalent to d(K bn , K kbn ).

The idea is then to find the largest b n , which minimizes this quantity. Several distances (Kolmogorov distance, Wasserstein metrics, etc...) may be used.

Of course, for large datasets, such a method is very computationally expensive. We recommend only choosing a limited range of values for b n and discretizing this range so as to compute the distance d(K bn , K kbn ) only on a limited number of points and to select the ones which minimize this quantity.

Other empirical approaches have been proposed to deal with the problem of the high volatility of subsampling distributions for too large subsampling sizes (see for instance [START_REF] Politis | Subsampling[END_REF] chap. 9). The idea is simply to look at the quantiles of subsampling distributions and to find the largest value such that the quantile remains stable. Some empirical arguments to understand the principle underlying this idea are given in [START_REF] Bertail | Comments on: Subsampling weakly dependent time series and application to extremes[END_REF]. Indeed a subsampling distribution may simply be seen as a U -statistic but with varying kernel of increasing size b n and its quantiles the inverse of a U-process. The main tools for studying the behavior of subsampling distribution are Hoeffding decomposition of the U -statistic and empirical process theory as considered in [START_REF] Arcones | Limit theorems for U-processes[END_REF]; [START_REF] Heilig | Limit theorems for the infinite-degree u-process[END_REF]. The difficulty in choosing the subsampling size is that, in comparison to U -statistics with a fixed degree, the linear part of the U -statistic is not always the dominating part in the Hoeffding decomposition. For rather small or moderate b n , it can be shown that the U -statistic is asymptotically normal with a convergence rate of order n bn . However, when b n becomes too large, the remainder in the Hoeffding decomposition dominates and the U -statistic behaves very erratically. This explains why the quantiles of subsampling distributions behave erratically as soon as b n is too large.

Subsampling in a growing environment

In many fields, data are now collected online. As a consequence, the size of the database may evolve quickly in time. Then, one may wish to update previous estimations without accessing the whole database again. How is it possible to update the subsampling distributions taking into account the new incoming data, when the size of the database is large and increases? To solve this problem, we present a very simple sequential algorithm.

The idea is as follows: assume that at time t, we have obtained a subsample without replacement of size b n (uniformly) from the original data. That is the probability of a given subsample is n bn -1 . At the time t + 1, the new sample size is n + 1. Then for this newcomer, proceed as follows:

• Draw a Bernoulli r.v. Z with parameter 1 -b n /(n + 1) • keep the original subsample if Z = 1, with probability 1 -b n /(n + 1),

• else with probability b n /(n + 1), choose one element of the current subsample (without replacement, uniformly with probability 1/b n ) and replace it with this newcomer.

If several newcomers arrive at the same time, then use sequentially the same algorithm by increasing the size of the population. Notice that this algorithm may be easily implemented sequentially to update all the subsamples already obtained by Monte Carlo simulations at some given time.

The arguments below show that the resulting algorithm is the realization of subsampling without replacement from the total new population.

It may be simply proved by recurrence [START_REF] Mcleod | A convenient algorithm for drawing a simple random sample[END_REF]. Indeed, assume that the probability of the original sample is n bn -1 then

• if Z = 1, the probability of the new sample is n bn -1 × (1 -bn n+1 ) = n+1 bn -1
• if Z = 0, the probability of the new sample is n bn

-1 × bn n+1 ( 1 bn + n-bn bn ) = n+1 bn -1 .
It follows that the corresponding subsample at any step is actually a subsample obtained without replacement from the total current population.

If we want to increase the size of the subsample, starting from a subsample of size b n in a population of size n then we simply draw uniformly without replacement from the n -b n remaining outside sample, observations with probability 1/(n -b n ). It may be sometimes easier (for instance using Apache Spark) to use sampling with replacement. It is known in that case that when b n is small enough such that bn √ n → 0, then the probability to draw the same individual twice converges to 0. Indeed, when bn √ n → 0, by Stirling formula, we have ( n bn ) n bn → 1, so that with and without replacement samplings are asymptotically equivalent under this condition.

Some empirical results

Simulation results

In this simulation section, the implementations were executed under R on a standard PC with a 5GHz Intel processor and 2G of Ram. The purpose is to show that we might gain a lot in terms of computation times for a large database and that using our confidence intervals for risks in a simple pattern recognition framework makes sense.

Maximum likelihood estimation for a simple logistic model

The purpose of this example is to explore the feasibility and the computer performance of the procedures described before in an estimation framework. We consider here a very simple toy model to highlight some inherent difficulties with subsampling. Consider a linear logistic regression model with parameter θ = (β 0 , β) ∈ R × R d . Let X be a d-dimensional marginal vector of the input random variables. The linear logistic regression model related to a pair (X, Y ) can be written as

P θ {Y = +1 | X} = exp(β 0 + β T X) 1 + exp(β 0 + β T X) .
In high-dimension i.e., when d is very large and for very large n, the computation of the full parametric maximum likelihood estimator (MLE) of θ may be difficult to obtain in a reasonable time. We assume that d << n but also that the subsampling sizes which will be used are such that d << b n .

For unbalanced populations (a lot of 1's in comparison with 0's and vice versa), the probability to get a subsample with only unit values (or zeros) is high and the MLE will not be convergent (a similar problem appears if the labels are fully separated). This is by no means contradictory with the asymptotic validity of subsampling in this case: it has been shown in Le [START_REF] Cam | Maximum likelihood: An introduction[END_REF] that the true variance of the MLE in a finite population is +∞. Subsampling simply reproduces this fact on a smaller scale.

In that case, one should condition on the fact that the ratio of the numbers of 1's by the number of 0's is not too small (or not too close to 1). Else, the subsample should be eliminated. We fix this ratio to 3% in our simulations.

In the following, we simulate the toy logistic model

Y i = 1 if 3X i + ε i > 0 0 else
with X i ∼ N (0, 1) and ε i independent logistic random variables with mean 0 and variance 1. We choose respectively n = 10 6 and n = 10 7 . We present the results for two randomly generated data set of large size, n (with one regressor variable here) and Table 1 shows the averages and sample variances of B = 999 subsample copies of the MLE (i.e., for fitting a logistic model) based on subsample draws of size b n and b n /e (when estimating the rate of convergence). In this example, the statistics T n to which we apply our results is thus the MLE of the parameter β Even for reasonable sizes, our estimation procedure proved to be useful. For instance in R, with 1 GB of memory, the usual libraries (sampleSelection, glm) fail to estimate the model with a size of n = 10 7 observations (for capacity reasons), whereas it takes only 12s to get subsampling based confidence bound with B = 999 replications of the procedure and b n = n 1/3 . Here, it is not required to estimate the rate of convergence since the rate τ n = n 1/2 is known, but we did it and present the result with an estimated rate with J = 2 subsampling distributions as proposed in this paper, one with size b n and the other with size b n /e. If we estimate the rate of convergence with J = 29 subsampling distributions based on subsampling sizes equal to n 1/3+j/(3(J-1)) , j = 0, . . . , 28, as done in the paper [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF], then one gets quite similar results in terms of variance of the MLE but with 999 × 29 simulations : it then takes around 14 times more time to complete these tasks on the same computer. We notice that the mean trick (pooled estimator) does not improve the MLE theoretically but allows a quicker and more feasible computation of the recentering factor.

The mean of the estimations of β (and the variances) over the B = 999 repetitions with the subsampling procedure are given in Table 1 for different values of b n n 1/3 , n 1/2 , n 2/3 given in the table and on the whole sample with the corresponding total execution times. v ar 1/2 gives a rescaled estimator of the variance (to compare with the estimator of the variance on the whole database). That is, we compute the variance of the subsampling distribution say s 2 bn,n and then extrapolate the variance to the original size by computing τ bn s 2 bn,n /τ n . When the rate is unknown the same trick is applied with an estimated rate. The estimated variance v ar(β n ) 1/2 ) is here obtained using the Hessian of the likelihood function at the MLE using the glm function: 10 7 is the order of the largest size for which we have been able to obtain estimator of this quantity in R with this simple model. Notice that even with a size of n 1/3 , we were able to get the correct order for the variance. In estimating β = 3, the bias in Table 1 can be large for small subsampling sizes but almost vanishes for n 2/3 . For this size, we get the same order as the one on the MLE on the whole database : but in terms of computation, n 2/3 is already too big, since in that case we were able to proceed with the MLE on the whole database in less than 5 minutes (whereas it takes 12 minutes to replicate 999 times the procedure on the n 2/3 sample size). But for n 1/2 , we get a gain of 4 (and 12 for n 1/3 ) for a similar accuracy: of course, this strongly depends on the degree of accuracy that one wishes to obtain on the parameter of interest and on the capacity of the computer.

Estimation of the out-of-sample error with k-nearest-neighbor algorithm.

Consider a simple pattern recognition framework. Assume that {(X 1 , Y 1 ), . . . , (X N , Y N )} is a sample of i.i.d. random pairs taking their values in some measurable product space X × {-1, +1}. In this standard binary classification framework, the multidimensional r.v.'s X are used to predict the binary label Y . The distribution P can also be described by the pair (F, η) where F (dx) denotes the marginal distribution of the input variable X and η(x) = P{Y = +1 | X = x}, x ∈ X , is the conditional distribution. The goal is to build a measurable classifier φ : X → {-1, +1} with minimum risk defined by

L(Y, φ(X)) def = I{φ(X) = Y }, (10)
where I{.} is the indicator function. It is well-known that the Bayes classifier φ * (x) = 2I{η(x) > 1/2} -1 is a solution of the risk minimization problem over the collection of all classifiers defined on the input space X .

It is now possible to apply the different subsampling procedures to different classes of functions (or algorithms) to estimate their prediction capability as described in part 3.3. We focus here on the out-of-sample errors. In this simulation part, we propose the use of F 1 : a parametric logit model, F 2 : the k-nearest neighbor method. We will also consider random forest models, SVM and neural networks in the case studies of part 4.2.

It is known that, under some regularity assumptions, the four methods are consistent so that, asymptotically, the approximation error ∆ n converges to 0 a.s.. For F 1 , it is known that the rate of convergence of the empirical risk is √ n; however for the other algorithms, even if some bounds exist on the generalization capability, the rate of convergence is not clear. Our method will allow us to evaluate the different algorithms and the rate of convergence of the algorithms.

Considering the same simulated data as in 4.1.1, we now use the subsampling method to estimate the out-of-sample errors of k-nearest neighbor (KNN) method and logit models on several subsampling sizes and compare them to that obtained on the full database. We consider a training set equal to 0.7n and a test set of size 0.3n (similar results have been obtained for other test sets). The computation times in Table 2 clearly show the computation gains. A striking result is for n = 10 7 because it takes almost 5 hours to get an estimator of this quantity on the whole sample whereas the subsampling method takes (at worst) 15 minutes with n 2/3 . Table 2 presents the average of out-of-sample errors (over B=999 subsamples) of the subsampling distribution. It seems that even with a size of order n 1/3 we still get a good approximation in less than 45 seconds. With the subsampling method by using an extrapolated variance as described in 4.1.1, we are also able to estimate the variance of the out-of-sample error (given in parenthesis in the table over the estimation on the whole sample). To evaluate the performance of our confidence intervals, we have repeated the preceding procedure M=500 times. First, for independent realizations (L=999), we evaluate the risk over the whole dataset and average this risk to have an idea of the true risk (average in Table 3). Then for m=1 to M=500 different samples, we construct confidence intervals for the risks of two methods: the one associated with the logistic regression and the one associated with KNN. We give in Table 3 the average of the lower (Lower) and upper (Upper) bounds of the confidence intervals. The coverage (cover) probability of our intervals is obtained by computing the number of times the "true risk" associated with the procedure belongs to the confidence intervals.

Table 3 shows that the intervals constructed with an interpolated variance using the average as re-centering factor is quite large in terms of errors even for quite large sample sizes and seem to be a little bit too conservative. However in this example we obtain narrower and clearly more accurate results for the Logit model which seems reasonable, since the data itself is generated in this true model. These simulations show that it is possible to compare in a reasonable time the outof-sample errors for several competing methods (with confidence intervals). We should mention that the computational time is not reasonable for neural networks because they require a long time to be trained and the parameters should be well-calibrated for each simulation. Despite this problem we also have implemented also this method in the case studies below.

Two case studies on real data sets

In this section, the implementation is performed in Python using the libraries NumPy, SciPy, Sklearn, Tensorflow, on an Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz machine with 144 cores and 250GB RAM.

Tested models

Subsampling techniques are implemented on potentially time-consuming procedures (Decision Tree DT [START_REF] Breiman | Classification and regression trees. Routledge Bretagnolle J (1983) Lois limites du bootstrap de certaines fonctionnelles[END_REF], Random Forest RF [START_REF] Breiman | Random forests[END_REF], Support Vector Machine SVM [START_REF] Chang | Libsvm: a library for support vector machines[END_REF], Neural nets (3 types: NeuralNet which is a fully connected multi-layer perceptron with one hidden layer [START_REF] Rumelhart | Learning representations by backpropagating errors[END_REF], NeuralNet3 which is a deeper multi-layer perceptron with three hidden layers, ConvNet (LeCun et al., 1989) which is a special architecture of a neural net that takes account of the hierarchical pattern in data and is commonly used in computer vision), a logit model Logit [START_REF] Mccullagh | Generalized linear models[END_REF], and K-nearest neighbors KNN.

Note that ConvNet is only used on MNIST data set as it is mainly used on image processing and that no hierarchical pattern exists between the features of VeReMi dataset.

The hyperparameters of all tested models were not specifically optimized to do the task on both data sets. We kept the default values found in sklearn that are based on recommendations from the original authors.

Description of the datasets

Vehicular Reference Misbehavior (VeReMi)

The VeReMi (Vehicular Reference Misbehavior) data set (see [START_REF] Van Der Heijden | VeReMi: A dataset for comparable evaluation of misbehavior detection in vanets[END_REF]) contains data about the detection of misbehavior in vehicular networks, a particularly sensitive topic in cooperative autonomous driving. The purpose of this application is to compare the relative risks of several algorithms and their confidence intervals. The VeReMi data set comprises N = 424, 810 individuals and only 12 variables.

EMNIST digits

The well-known EMNIST data set studied for instance in [START_REF] Lecun | Comparison of learning algorithms for handwritten digit recognition[END_REF] contains binary images of handwritten digits and the corresponding true digits; the wellknown purpose is to find an algorithm that correctly recognizes the digit, see [START_REF] Cohen | EMNIST: an extension of MNIST to handwritten letters[END_REF]. The EMNIST data set comprises N = 240, 000 images and 784 variables (28 × 28 pixels for each image).

A smaller version of this database is the MNIST digits popularized by Lecun. The data set comprises N = 60, 000 images. This version of the database is used to check whether the estimation of the convergence rate is similar when based on a smaller data set. We refer to this database as "Lecun MNIST data set".

Comparison of the performances of the tested models

Computing Out-of-sample errors Figures 1, 2 and 3 provide the estimation of the out-of-sample errors of several models with 90% confidence intervals as well as the computing times for each of them, according to subsampling sizes, for the VeReMi and EMNIST digits data sets. The methodology is described in sections 3.1 and 3.2: B is set to 999, b n is ranging from N 1/3 to N 2/3 similarly to what we did in the previous section, with the estimation of J = 29 subsampling distributions. For VeReMi, we observe the superiority of the two tree-based approaches (RF and DT), and the lower performance of the KNN approach in terms of errors, but DT would clearly be preferred as its computing time does not explode as that of RF.

For the EMNIST digit, most confidence intervals overlap for subsampling sizes up to b n = 4, 000 although ConvNet and SVM show the lowest errors. From this error plot, we can not conclude that their superiority is significant. In terms of computing times, SVM would clearly be preferred to ConvNet as the computing time for the latter is greater than 6 hours for b n = 4, 000.

For Lecun MNIST data set, only 6 models were compared, and SVM and RF have the best performances in terms of errors, with SVM showing lower computing times. Note that the computing times on this data set 4 times smaller are also 4 times better.

The RF, NeuralNet and SVM models have a higher execution time than other models for all subsampling sizes. In particular, for the VeReMi data set where the number of variables is very small, NeuralNet has a relatively smaller execution time compared to Random Forests due to the reduced number of parameters in the network. This changes for the EMNIST digit, where the number of features is bigger and Neural-Net takes longer to train. The simpler models KNN, DT, and Logit take relatively smaller execution times. Overall all execution times increase with subsampling sizes with different rates for each model. Figures 4 shows the extrapolation to the full dataset sizes of the out-of-sample errors. It requires the computation of the estimated convergence rate presented in the next section, that is, the computation of the J = 29 subsampling distributions (the results with only J = 2 are quite similar).

For VeReMi, RF, and DT perform the best, similarly to each other. For EMNIST, Convnet clearly outperforms the other models, SVM is next. For Lecun MNIST, SVM, and RF do not differ much and we can see that NeuralNet behaves better on the N = 240, 000 EMNIST data set than the 60, 000 Lecun MNIST data set.

In comparison with the Out-of-sample errors on the full dataset (Table 4, 5), the order of best to worst performing methods matches the one with the extrapolated error. Except for KNN on the VeREMi data set where the introduction of the full data helped achieve much better performance. We estimate here the rate of convergence τ n as n α using 2 and 29 different values of subsampling size: in terms of estimation, the rate is of the same order. See figures 5 to 11 to understand why: just pick two points (not too close) among all subsampling sizes and it somehow gives the right slope. However, it seems that using more subsampling distributions as in [START_REF] Bertail | Subsampling the distribution of diverging statistics with applications to finance[END_REF] gives more precise results (we thus only give these estimators in Table 6 and7). As described at the end of section 2.3, we consider here a regression of log range, on the log of the subsampling size. We propose 3 different ranges, either the inter-quartile-range (percentiles at 25% and 75%) denoted IQR, or an inter-percentile range with length 80% denoted IPR80 (percentiles at 10% and 90%), or with length 90% denoted IPR90 (percentiles at 5% and 95%). Individual graphs showing the slope α of each model/IPR choice are postponed to the appendix (see Figures B1,B2,B3,B4,B5,B6,B7) and results are summarized in Tables 6 and7. For both data sets, we observe that the results are quite stable when changing the IPR approach as can be seen from the figures in the appendix, or using only 2 subsampling distributions. For VeReMi, we observe that RF and DT have faster convergence rates (α ≈-.65) than the 4 other models (α ≈ -1/2). These models also reached the lowest Out-of-sample errors. For EMNIST, we tested two more models: a Neural Net with 3 layers and a Convolution Net to see how the rates of convergence would be impacted. We observe that here, Logit, DT and NeuralNet3 models have α ≈ -1/2 while KNN and RF have α ≈ -2/3, and SVM and ConvNet have α ≈ -3/4, NeuralNet being closer to α ≈ -.6. Again, the models with faster convergence rates are also those with lower out-of-sample errors. However, we see that the convergence rates do depend on the learning problem as we observe for instance that DT is very efficient in the VeReMi case (12 features) and not in the EMNIST case (784 features). . Now for the last propositions, just notice that we have the decomposition K bn (x | X n , τ . ) -K(x, P ) = K bn (x | X n , τ . ) -K bn (x, P ) + K bn (x, P ) -K(x, P ) and K bn (x | X n , τ . ) -K n (x, P ) = K bn (x | X n , τ . ) -K(x, P ) + K(x, P ) -K n (x, P ) Now use assumption A5 and the preceding results to conclude.

Estimation of the convergence rates

A.2. Proof of Lemma 1 from section 2.3

For any > 0, we know from Theorem 1 in section 2.2, that there exists some L = L ε , Pr P {|K bn (x | X n , τ . ) -K(x, P )| ≥ L/δ β (n)} ≤ ε (A1) uniformly in x. Put η n = L δβ(n) and define the quantile z = K -1 bn (t -η n | X n , τ . ), then K bn (z | X n , τ . ) ≥ t -η n . Combining this with (A1) implies that z ≥ K -1 (t -2η n , P ).with probability at most ε. Similarly, define y = K -1 (t, P ), then we have y ≥ K -1 bn (t -η n | X n , τ . ) with probability at most ε. Hence, for any t and any > 0, we have the inequality :

Pr P K -1 (t -2η n , P ) ≤ K -1 bn (t -η n | X n , τ . ) ≤ K -1 (t, P ) ≤ 2ε. ( A2
)
This clearly yields , for → 0 + , that K -1 bn (t | X n , τ . ) = K -1 (t, P ) + o P (1).

But, by assumption A6 and using the Inverse function theorem for Lipschitz, strictly increasing continuous functions that, there exists an L ε such that |K -1 (t -2η n , P ) -K -1 (t, P )| ≤ L ε η n

Using again (A2) twice, we get that K -1 bn (t | X n , τ . ) = K -1 (t, P ) + O P (η n ) = K -1 (t, P ) + O P (δ β (n) -1 ), uniformly in the neighborhood of t.

A.3. Proof of Theorem 2 from section 3.1

The proof follows the same lines as [START_REF] Bertail | On subsampling estimators with unknown rate of convergence[END_REF]. For any x, consider

K bn (x | X n , τ . ) ≡ q -1 q i=1 1{b α n (T bn,i -θ n ) ≤ x} = q -1 q i=1 1{b α n (T bn,i -θ) -b α n ( θ n -θ) ≤ x}
and define the correctly recentered U -statistic Conditions A2 imply that Pr(E n ) -→ n→∞ 1; hence, with probability tending to one, we get that

U n (x) = q -1 q i=1 1{b α n (T bn,i -θ) ≤ x}
U n (x -) ≤ K bn (x | X n , τ . ) ≤ U n (x + ).
Let us show that U n (x) converges to K(x, P ) in probability. For this, introduce the U -statistic with varying kernel defined by : V n (x) = q -1 q i=1 1{b α n (T bn,i -θ) ≤ x}, which is the equivalent of U n (x), with the true rate rather than the estimated one.

Recall that since V n (x) is a U -statistic of degree b n , such that bn n → 0, by Hoeffding's inequality, we have V n (x) = K(x, P ) + O P (1/ (b n /n)) as n → ∞ in probability. Now, for any 1 > 0, we have that, with probability tending to 1,

U n (x) = q -1 q i=1 1 b α n (T bn,i -θ) ≤ b α n b α n
x ≤ V n (x + 1 ).

A similar argument shows that we also have U n (x) ≥ V n (x -1 ) with probability tending to one. But we have V n (x + 1 ) → K(x + 1 , P ) and V n (x -1 ) → K(x -1 , P ) in probability. Therefore, letting 1 → 0, we have that U n (x) → K(x, P ) in probability as required.

Proving that we have K n (x, P ) -K(x, P ) → 0 as n → ∞.

follows now by the same arguments as before by recalling that K n (x, P ) = Pr τ n (T n -θ) ≤ x τ n τ n = Pr(τ n (T n -θ) ≤ x(1 + o P (1)) and using the continuity of the limiting distribution.

The second part of the theorem is a straightforward consequence of the uniform convergence of K bn (x | X n , τ . )-K n (x, P ) to 0, over the neighborhood of any continuity point of the true limiting distribution.

The last result is obtained by the same arguments, by just replacing by log(n)δ β (n) -1 . In that case, the probability of the event E n is controlled by assumption A4. All the approximations o P (1) then becomes O P (log(n)δ β (n) -1 ) using the same arguments as in Theorem 1. Similarly 1 can be replaced by 1 log(b n )δ β (n) -1 which is anyway smaller than O P (log(n)δ β (n) -1 ).

Proof of corollary 1 from section 3.3.1

Recall that θ q n = 1 q 1≤j≤q E (j)

bn . Notice first that the value E bn | < ε}. Then, by the preceding Hoeffding's inequality, using the fact that B = n γ , we get :

Pr(E c bn ) ≤ exp - 2ε 2 (n -b n ) τ 2 bn M 2 + γ ln(n) , (A3)
which goes to 0 under our assumptions on b n . It follows that Pr(E bn ) → 1 as n → ∞.

As in the proof of Theorem 3, we have a Bernstein inequality for τ bn | θ q n -θ * | . Now, apply the same arguments as in Theorem 3, with T bn,i = E which is exactly the one of the U -statistics of Theorem 1, so that the preceding results apply.

  divide as well this set into a training set D (k) bn,T r of size p 0 b n and a test set D (k) bn,T est of size (1 -p 0 )b n to construct the corresponding estimator

Figure 1 .Figure 2 .

 12 Figure 1. Comparison of Out-of-sample errors and their associated computing times from 6 different models according to the subsampling size, VeReMi data set. For errors, 90% confidence intervals are provided.

  and the eventE n = {b α n | θ n -θ| ≤ },for some > 0.Since α = α+o P ((log n) -1 ), we have as well n α = n α (1+o P (1)) and b α n = b α n (1+o P (1)). Notice for the last statement of the theorem, that ifα = α + O P (δ β (n) -1 ), we get n α = n α (1 + O P (log(n)/δ β (n))) and b α n = b α n (1 + O P (log(b n )/δ β (n)))).

  n that can be controlled by standard arguments on sums. Indeed, by Hoeffding's inequality, we have that, for some constant M > 0,bn | > x ≤ 2B exp -2x 2 (n -b n ) M 2 .Now, notice that the subsampling distribution may be writtenK (B) bn (x | X n , τ . ) = B -1 bn + θ * -θ q n .As in the proof of Theorem 1, consider the event E bn = {τ bn sup

.

  bn to get that, some constants C1, C 2 and for our b n , for any fixed ε > 0Pr τ bn n b n | θ q n -θ * | > ε ≤ C 1 exp -C 2 ε 2 + Pr(E c bn )From this and equation (A3), we get that| θ q n -θ * | = O P τ -1 bn bn nThis result and again equation (A3) imply that the τ bn a (j)n 's are uniformly small. Using the continuity of the limiting distribution as in the proof of Theorem 1, it follows that it is sufficient to study the subsampling distribution K bn -θ * ) ≤ x}.

Table 1 .

 1 MLE variance estimations for a logistic model with n = 10 6 , 10 7 .

		subsample (B = 999 replications)	whole sample
						βn
	n	b n	βbn	v ar 1/2	time	( v ar(β n ) 1/2 )
						time
		n 1/3 ≈ 100	3.19 0.0064	13 s	2.992
	10 6	n 1/2 ≈ 1000 3.022 0.0063	36 s	(0.0061)
		n 2/3 ≈ 10000 2.996 0.0060 3.26 mn	28.75 s
		n 1/3 ≈ 215	3.10 0.0020	41 s	2.998
	10 7	n 1/2 ≈ 3162 3.009 0.0020 1.25 mn	(0.0019)
		n 2/3 ≈ 46415 2.998 0.0019	12 mn	4.69 mn

Table 2 .

 2 Estimation of the out-of-sample error by subsampling and on the whole sample -KNN model

	KNN	subsample (B = 999 replications)	whole sample
	n	b n	out-of-samp. error	time	out-of-samp. err	time
		n 1/3	0.1177	4.79 s	0.1158	5.252 mn
	10 6	n 1/2	0.1165	5.76 s	(0.008)
		n 2/3	0.1167	43.5 s
		n 1/3	0.1166	44.7 s	0.1141	4h57mn
	10 7	n 1/2	0.1163	50.7 s	(0.006)
		n 2/3	0.1161	15.35 mn

Table 3 .

 3 Risk, confidence intervals and coverage probability: KNN and logit models

	M ethod KN N				Logit
	size	Lower U pper average cover Lower U pper average cover
	n = 10 6 0.032	0.207	0.118	0.98	0.012	0.110	0.066	0.96
	n = 10 7 0.046	0.185	0.115	0.96	0.017	0.108	0.061	0.96

Table 4 .

 4 Summary statistics of Out-of-sample errors from different models on the full EMNIST data set size

		Logit	RF	SVM	DT	KNN NeuralNet
	count	10	10	10	10	10	10
	mean 0.062435 0.019450 0.011373 0.080225 0.018935	0.028037
	std	0.000742 0.000503 0.000431 0.001078 0.000345	0.000926
	min	0.061417 0.019000 0.010604 0.078000 0.018542	0.026792
	10%	0.061623 0.019056 0.010792 0.078844 0.018617	0.027110
	50%	0.062375 0.019208 0.011385 0.080719 0.018844	0.028010
	90%	0.063517 0.020296 0.011827 0.081125 0.019500	0.028950
	max	0.063667 0.020333 0.011958 0.081125 0.019500	0.029812

Table 5 .

 5 Summary statistics of Out-of-sample errors from different models on the full VeReMi dataset size

		logit	RF	SVM	DT	Knn NeuralNet
	count	10	10	10	10	10	10
	mean 0.254134 0.000194 0.254134 0.001743 0.030807	0.254228
	std	0.001391 0.000050 0.001391 0.000208 0.000526	0.001340
	min	0.251442 0.000141 0.251442 0.001565 0.029849	0.251677
	10%	0.252957 0.000152 0.252957 0.001587 0.029933	0.252948
	50%	0.253878 0.000177 0.253878 0.001671 0.030920	0.254172
	90%	0.255643 0.000262 0.255643 0.002120 0.031286	0.255632
	max	0.255738 0.000294 0.255738 0.002130 0.031402	0.255738

Table 6 .

 6 Estimation of the convergence rates of the 6 different models, with 3 different methods -the VeReMi data set

	Model	IQR IPR80 IPR90
	NeuralNet -0.487 -0.500 -0.499
	Logit	-0.493 -0.513 -0.498
	KNN	-0.497 -0.506 -0.493
	SVM	-0.510 -0.494 -0.504
	RF	-0.595 -0.636 -0.635
	DT	-0.611 -0.624 -0.630

Table 7 .

 7 Estimation of the convergence rates of the 8 different models, with 3 different methods -the EMNIST data setTable 8 also compares the EMNIST IQR results with the Lecun MNIST IQR results. The smaller dataset shows somehow lower rates than the larger one.

	Model	IQR IPR80 IPR90
	Logit	-0.558 -0.560 -0.562
	NeuralNet3 -0.560 -0.514 -0.525
	DT	-0.567 -0.542 -0.533
	NeuralNet	-0.591 -0.587 -0.578
	KNN	-0.657 -0.655 -0.673
	RF	-0.688 -0.675 -0.687
	ConvNet	-0.722 -0.753 -0.750
	SVM	-0.773 -0.754 -0.757

Table 8 .

 8 Estimation of the convergence rates of the 6 different models (IQR method), comparison of EMNIST data set (N = 240, 000) and Lecun MNIST data set (N = 60, 000) and M 2 , we also have an exponential inequality for K bn (x | X n , τ . ) of the form Pr (|K bn (x | X n , τ . ) -K bn (x, P )| > ε) ≤ M 1 exp(-

	Model	EMNIST Lecun MNIST
	Logit	-0.558	-0.564
	RF	-0.688	-0.656
	DT	-0.567	-0.518
	KNN	-0.657	-0.641
	NeuralNet	-0.591	-0.427
	SVM	-0.773	-0.682

Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 Slope: 0.4999 -3.0 -2.5 -2.0 -1.5

Subsampling size (log) Log Inter Percentile Range 0.1, 0.9 Generalization error -NeuralNet Figure B2. Estimation of convergence rate, based on IPR80, VeReMi dataset, N = 424, 810

Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499 Slope: 0.499

Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869 Slope: 0.6869

Generalization error -SVM n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 B.1.3. IPR90 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -KNN n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -LOGIT n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -SVM n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 B.2. Detailed results on EMNIST digit dataset B.2.1. IQR n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -KNN n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -LOGIT n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -SVM n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -KNN n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 B.2.3. IPR90 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -KNN n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -LOGIT n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -RF n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 Generalization error -SVM n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 B.3. Detailed results on Lecun MNIST digit dataset B.3.1. IQR n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29 n_b: 29

This appendix is dedicated to the proof of the theorems and the corollary of this article.

A.1. Proof of Theorem 1 from section 2.2

Introduce the U -statistic

1{τ bn (T bn,i -θ) ≤ x}.

Then, we have the simple decomposition

Since E P [V bn (x)] = K bn (x, P ) and V bn is a U -statistic of degree b n with kernel bounded by 1, we have by Hoeffding's inequality

Now, we can write using the same argument (twice), for any η > 0, But since K bn (x, P ) is supposed to be continuous at x (at least asymptotically), for n large enough, the last term is 0 for a well chosen η. More precisely, under A5 and the Lipschitz condition A6, we have It follows that if we choose η such that η < ε/12L, for n large enough the last term vanishes. Now for this choice, we get that, by hypothesis (A4), for some non-negative