François Bachoc

Asymptotic analysis of covariance parameter estimation for Gaussian processes in the misspecified case

In parametric estimation of covariance function of Gaussian processes, it is often the case that the true covariance function does not belong to the parametric set used for estimation. This situation is called the misspecified case. In this case, it has been observed that, for irregular spatial sampling of observation points, Cross Validation can yield smaller prediction errors than Maximum Likelihood. Motivated by this comparison, we provide a general asymptotic analysis of the misspecified case, for independent observation points with uniform distribution. We prove that the Maximum Likelihood estimator asymptotically minimizes a Kullback-Leibler divergence, within the misspecified parametric set, while Cross Validation asymptotically minimizes the integrated square prediction error. In a Monte Carlo simulation, we show that the covariance parameters estimated by Maximum Likelihood and Cross Validation, and the corresponding Kullback-Leibler divergences and integrated square prediction errors, can be strongly contrasting. On a more technical level, we provide new increasing-domain asymptotic results for the situation where the eigenvalues of the covariance matrices involved are not upper bounded.

) also suggest that CV-type methods can provide better predictions than ML-type methods under misspecification.

In this paper, we aim at providing a general asymptotic analysis of the misspecified case, which would confirm the findings of the aforementioned references on the comparison between ML and CV. The two most studied asymptotic frameworks are the increasing-domain and fixed-domain asymptotics (Stein, 1999, p.62). In increasing-domain asymptotics, the average density of observation points is bounded, so that the infinite sequence of observation points is unbounded. In fixed-domain asymptotics, this sequence is dense in a bounded domain.

. Microergodic covariance parameters have an asymptotic influence on predictions, as shown in (Vazquez, 2005, Ch.5).

Nevertheless, similarly to Bachoc (2014), we do not address fixed-domain asymptotics in this paper. In Bachoc (2014), the reason is that fixed-domain asymptotics yields no impact of the spatial sampling of observation points on covariance function estimation. We follow this line here, since Bachoc (2013) makes the important observation that the comparison between ML and CV strongly depends on the spatial sampling. In this paper, we indeed derive asymptotic results for CV for which the type of spatial sampling used is crucial. Furthermore, another downside of fixed-domain asymptotics is that the results currently under reach, despite their significant insights, are restricted in terms of covariance model. For example, Ying (1993) addresses ML for the tensorized exponential model only and Loh (2005) addresses ML for the Matérn 3/2 covariance model only.

. This is because increasing-domain asymptotics is characterized by a vanishing dependence between observations from distant observation points, so that a large sample size gives more and more information about the covariance structure.

The increasing-domain asymptotic spatial sampling we consider consists of n independent observation points with uniform distribution on [0, n 1/d ] d , for d ∈ N * . We prove that ML asymptotically minimizes, within the misspecified model, the Kullback-Leibler divergence from the true covariance function, defined at the observation vector. As we discuss after Theorem 3.3, this can actually be shown under more general spatial samplings than that considered here, but still necessitates an original proof. In the misspecified case, the important point is that no information is provided on the quality of the ML estimator for subsequent predictions of the Gaussian process at new points.

, that when the spatial sampling is not too regular, CV generally performs better than ML for prediction in the misspecified case, while it can instead perform poorly under regular sampling. On a more

Introduction

Kriging models [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] consist in inferring the values of a Gaussian random field given observations at a finite set of observation points. They have become a popular method for a large range of applications, such as numerical code approximation [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] and calibration [START_REF] Paulo | Calibration of computer models with multivariate output[END_REF] or global optimization [START_REF] Jones | Efficient global optimization of expensive black box functions[END_REF].

One of the main issues regarding Kriging is the choice of the covariance function for the Gaussian process. Indeed, a Kriging model yields an unbiased predictor with minimal variance and a correct predictive variance only if the correct covariance function is used. The most common practice is to statistically estimate the covariance function, from a set of observations of the Gaussian process, and to plug (Stein, 1999, Ch.6.8) the estimate in the Kriging equations. Usually, it is assumed that the covariance function belongs to a given parametric family (see [START_REF] Abrahamsen | A review of Gaussian random fields and correlation functions[END_REF] for a review of classical families). In this case, the estimation boils down to estimating the corresponding covariance parameters.

For covariance parameter estimation, Maximum Likelihood (ML) is the most studied and used method, while Cross Validation (CV) [START_REF] Sundararajan | Predictive approaches for choosing hyperparameters in Gaussian processes[END_REF][START_REF] Zhang | Kriging and cross validation for massive spatial data[END_REF] is an alternative technique. The two estimators have been compared by several references. Consider first the case where the true covariance function of the Gaussian process belongs to the parametric family of covariance functions used for estimation, that we call the well-specified case. Then [START_REF] Stein | A comparison of generalized cross validation and modified maximum likelihood for estimating the parameters of a stochastic process[END_REF] shows that for the estimation of a signal-to-noise ratio parameter of a Brownian motion, CV has twice the asymptotic variance of ML. In the situations treated by [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF], the asymptotic variance is also larger for CV than for ML. Several numerical results, showing an advantage for ML over CV as well, are available, coming either from Monte Carlo studies as in (Santner et al., 2003, Ch.3) or deterministic studies as in [START_REF] Martin | On the use of Kriging models to approximate deterministic computer models[END_REF]. The settings of both the above studies can arguably be classified in the well-specified case, since the interpolated functions are smooth, and the covariance structures are adapted, being Gaussian in [START_REF] Martin | On the use of Kriging models to approximate deterministic computer models[END_REF] and having a free smoothness parameter in [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Finally, in situations similar to the well-specified case, ML-type methods have been shown to be preferable over CV-type methods in [START_REF] Stein | Spline smoothing with an estimated order parameter[END_REF] for estimation and prediction. technical level, the proof of Theorem 3.4 tackles a case where the eigenvalues of the covariance matrices involved are not upper bounded as n → ∞. To the best of our knowledge, this situation has never been addressed in the increasing-domain asymptotic literature.

We conclude this paper by a Monte Carlo simulation, illustrating Theorems 3.3 and 3.4. The simulation highlights that the ML and CV estimators can estimate radically different covariance parameters, and that their subsequent performances for the Kullback-Leibler divergence and the integrated square prediction error can be strongly contrasting.

The rest of the paper is organized as follows. We present the context on parametric covariance function estimation in the misspecified case and on the spatial sampling in Section 2. We give the asymptotic optimality results for ML and CV in Section 3. We discuss the simulation results in Section 4. All the proofs are given in the appendix.

Context

Presentation and notation for the covariance model

We consider a stationary Gaussian process Y on R d with zero mean function and covariance function K 0 . Noisy observations of Y are obtained at the random points X 1 , ..., X n ∈ R d , for n ∈ N * . That is, for i = 1, ..., n, we observe y i = Y (X i ) + i , where = ( 1 , ..., n ) t , Y and (X 1 , ..., X n ) are mutually independent and follows a N (0, δ 0 I n ) distribution, with δ 0 ≥ 0 and I n the identity matrix of size n. The distribution of (X 1 , ..., X n ) is specified and discussed in Condition 2.3 below.

The case where Y is observed exactly is treated by this framework by letting δ 0 = 0. Otherwise, letting δ 0 > 0 can correspond for instance to measure errors [START_REF] Bachoc | Calibration and improved prediction of computer models by universal Kriging[END_REF] or to Monte Carlo computer experiments [START_REF] Le Gratiet | Asymptotic analysis of the learning curve for gaussian process regression[END_REF]. Note also that the case of a Gaussian process with discontinuous covariance function at 0 (nugget effect) is mathematically equivalent to this framework if the observation points X 1 , ..., X n are two by two distinct. [This is the case in this paper, in an almost sure sense, see Condition 2. 3.] Let p ∈ N * and let Θ be the compact subset [θ inf , θ sup ] p with -∞ < θ inf < θ sup < +∞. We consider a parametric model attempting to approximate the covariance function K 0 and the noise variance δ 0 : {(K θ , δ θ ), θ ∈ Θ}, with K θ a stationary covariance function and δ θ > 0. We call the case where there exists θ 0 ∈ Θ so that (K 0 , δ 0 ) = (K θ0 , δ θ0 ) the well-specified case. The converse case, where (K 0 , δ 0 ) = (K θ , δ θ ) for all θ ∈ Θ is called the misspecified case.

The well-specified case has been extensively studied in the Gaussian process literature, see the references given in Section 1. Nevertheless, the misspecified case can occur in many practical applications. Indeed, even if we assume δ θ = δ 0 for all θ, the standard covariance models {K θ , θ ∈ Θ} are often driven by a limited number of parameters and thus restricted in some ways. For instance, one common practice [START_REF] Martin | On the use of Kriging models to approximate deterministic computer models[END_REF]) is to use the Gaussian covariance model, where p = d + 1, Θ ⊂ (0, ∞) p , θ = (σ 2 , 1 , ..., d ) and K θ (t) = exp(-

d i=1 t 2 i / 2 i ).
With the Gaussian covariance model, all the covariance functions K θ generate Gaussian process realizations that are almost surely infinitely differentiable. Thus, the Gaussian model is de facto misspecified if the realizations of Y have only a finite order of differentiability. In theory, the Matérn model considered in Section 4 brings a solution to this problem, by incorporating a tunable smoothness parameter ν > 0. However, it is also common practice to enforce a priori this parameter ν to a fixed value (e.g. 3/2 in [START_REF] Chevalier | Fast parallel krigingbased stepwise uncertainty reduction with application to the identification of an excursion set[END_REF]).

In this paper, we are primarily interested in analyzing the misspecified case and stressing that the conclusions for it differ from those that are usually derived from the well-specified case. Nevertheless, the asymptotic results that are given in Section 3 actually do not assume one of the two cases and are valid for both.

We let X = (X 1 , ..., X n ) be the random n-tuple of the n observation points. For θ ∈ Θ, we define the n × n random matrix R θ by (R θ ) i,j = K θ (X i -X j ) + δ θ 1 i=j . We define the n × n random matrix R 0 by (R 0 ) i,j = K 0 (X i -X j ) + δ 0 1 i=j . We do not write explicitly the dependence of R θ and R 0 with respect to X and n. We define the random vector y = (y 1 , ..., y n ) t of size n by y i = Y (X i ) + i . Then, conditionally to X, y follows a N (0, R 0 ) distribution and is assumed to follow a N (0, R θ ) distribution under parameter θ. We do not write explicitly the dependence of y with respect to X, n and .

Maximum Likelihood and Cross Validation estimators

The Maximum Likelihood (ML) estimator is defined by θML ∈ argmin θ L θ , with

L θ := 1 n log (det (R θ )) + 1 n y t R -1 θ y. (1)
Remark 2.1. For concision, we do not write explicitly the dependence of L θ on X, n, Y and . We make the same remark for the CV criterion in (2) and (3).

Remark 2.2. In this paper, we enable the criterion (1) to have more than one global minimizer, in which case, the asymptotic results of Section 3 hold for any sequence of random variables θML minimizing it. The same remark can be made for the CV criterion (2). We refer to Remark 2.1 in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] for the existence of measurable minimizers of the ML and CV criteria.

The functional L θ is the modified opposite log-likelihood. It is globally admitted that, in the well-specified case, the ML estimator is preferable over most other potential estimators, both from asymptotic grounds and in practice. Under increasing-domain asymptotics, ML is consistent and asymptotically normal, with mean vector 0 and covariance matrix the inverse of the Fisher information matrix. This is shown in [START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF], assuming either some convergence conditions on the covariance matrices and their derivatives or gridded observation points. Similar results are provided for Restricted Maximum Likelihood in Cressie andLahiri (1993, 1996). In [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] asymptotic normality is also shown for Maximum Likelihood, using only simple conditions on the covariance model and for observation points that constitute a randomly perturbed regular grid.

The Cross Validation (CV) estimator, minimizing the Leave One Out (LOO) mean square error is defined by θCV ∈ argmin θ CV θ , with

CV θ := 1 n n i=1 (y i -ŷi,θ ) 2 , (2)
where ŷi,θ := E θ|X (y i |y 1 , ..., y i-1 , y i+1 , ..., y n ) is the LOO prediction of y i with parameter θ. The conditional mean value E θ|X denotes the expectation with respect to the distribution of Y and with the covariance function K θ and the variance δ θ , given X.

Let r i,θ = (K θ (X i , X 1 ), ..., K θ (X i , X i-1 ), K θ (X i , X i+1 ), ..., K θ (X i , X n ))
t . Define r i,0 similarly with K 0 . Define the (n -1) × (n -1) covariance matrix R i,θ as the matrix extracted from R θ by deleting its line and column i. Define R i,0 similarly with R 0 . Then, with y -i = (y 1 , ..., y i-1 , y i+1 , ..., y n ) t , we have ŷi,θ = r t i,θ R -1 i,θ y -i . Note that ŷi,θ is invariant if K θ and δ θ are multiplied by a common positive constant. Thus, the CV criterion (2) is designed to select only a correlation function K θ /K θ (0) and a corresponding relative noise variance δ θ /K θ (0). For the question of selecting the variance K θ (0) by CV, two cases are possible, depending on the parametric model:

First, several parametric models can be written so that they satisfy the decomposition θ = σ 2 , θ , with

σ 2 > 0, θ ∈ R p-1 and (K θ , δ θ ) = (σ 2 Kθ , σ 2 δθ )
, with Kθ a stationary correlation function and δθ ≥ 0. Hence, in this case, θ would be estimated by (2) in a first step. In a second step, σ 2 can be estimated as in [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF] by the equation σ2

CV ( θ) = (1/n) n i=1 {y i -ŷi, θ } 2 /c 2 i, θ
, where c 2 i, θ = var θ|X (y i |y 1 , ..., y i-1 , y i+1 , ..., y n ) is the LOO predictive variance for y i with parameters σ 2 = 1 and θ. The conditional variance var θ|X denotes the variance with respect to the distribution of Y and with the covariance function Kθ and the noise variance δθ , given X. Even if the decomposition θ = σ 2 , θ is possible, we address only the estimation of θ in this paper, that is given by (2).

Second, assume that the the aforementioned decomposition is not possible. One important example is the one addressed in Section 4, where a fixed value δ 1 of the noise variance is assumed, so that δ θ = δ 1 for all θ ∈ Θ. Then, usually, there is a unique parameter θCV minimizing (2). In this case, the estimated (K θCV , δ θCV ) is relevant in the misspecified case only for providing the predictor E θCV |X (Y (t)|y) of the value of Y at a new point t. Other quantities like var θCV |X (Y (t)|y), where var θ|X denotes the variance under parameter θ given X, can be unreliable. Nevertheless, the predictor E θCV |X (Y (t)|y) alone has the same applicability as many regression techniques like kernel regression or neural network methods and can be used in a wide range of practical applications. Though one can consider other Cross Validation criteria to optimize, like the log predictive probability [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], chapter 5, [START_REF] Zhang | Kriging and cross validation for massive spatial data[END_REF], [START_REF] Sundararajan | Predictive approaches for choosing hyperparameters in Gaussian processes[END_REF]), this paper is dedicated to the criterion (2), and toward analyzing it according to the accuracy of the predictor E θCV |X (Y (t)|y).

The two aforementioned cases are treated in the same fashion in the rest of the paper. Indeed, we only work under the assumption that θCV is a random variable minimizing (2) which is allowed to have more than one global minimizer (Remark 2.2).

The criterion (2) can be computed with a single matrix inversion, by means of virtual LOO formulas (see e.g [START_REF] Ripley | Spatial Statistics[END_REF]; [START_REF] Dubrule | Cross validation of Kriging in a unique neighborhood[END_REF]). These virtual LOO formulas yield, when writing diag (A) for the matrix obtained by setting to 0 all the off diagonal elements of a square matrix A,

CV θ := 1 n y t R -1 θ diag R -1 θ -2 R -1 θ y, (3)
which is useful both in practice (to compute the CV criterion quickly) and in the proofs for CV. Finally, in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] it is shown that, in the well-specified case, the CV estimator is consistent and asymptotically normal for estimating correlation parameters, under increasing-domain asymptotics. Nevertheless, in the well-specified case, ML is preferable as discussed above. In the misspecified case, this hierarchy can change, as shown numerically in [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF]. In Section 3, we give asymptotic grounds to this finding.

Random spatial sampling

We consider an increasing-domain asymptotic framework where the observation points are independent and uniformly distributed.

Condition 2.3. For all n ∈ N * , the observation points X 1 , ..., X n are random and follow independently the uniform distribution on [0, n 1/d ] d . The three variables Y , (X 1 , ..., X n ) and are mutually independent. Condition 2.3 constitutes an increasing-domain asymptotic framework in the sense that the volume of the observation domain is n and the average density of observation points is constant. Some authors define increasingdomain asymptotics by the condition that the minimum distance between two different observation points is bounded away from zero (e.g. [START_REF] Zhang | Toward reconciling two asymptotic frameworks in spatial statistics[END_REF]), which is not the case here. In [START_REF] Lahiri | Central limit theorems for weighted sums of a spatial process under a class of stochastic and fixed designs[END_REF] and [START_REF] Lahiri | Asymptotic distributions of M-estimators in a spatial regression model under some fixed and stochastic spatial sampling designs[END_REF], the term increasing-domain is also used, when points are sampled randomly on a domain with volume proportional to n. Observation points that are sampled uniformly and independently constitute the archetype of an irregular spatial sampling. Thus, as discussed in Section 1, CV is expected to perform well under Condition 2.3. This is indeed the case, as is shown in Section 3.

Asymptotic optimality results

Technical assumptions

We shall assume the following condition for the covariance function K 0 , which is satisfied in all the most classical cases, and especially for the Matérn covariance function. Let |t| = max i=1,...,d |t i |.

Condition 3.1. The covariance function K 0 is stationary and continuous on R d . There exists

C 0 < +∞ so that for t ∈ R d , |K 0 (t)| ≤ C 0 1 + |t| d+1 .
Next, the following condition for the parametric set of covariance functions and noise variances is slightly non-standard but not restrictive. We discuss it below.

Condition 3.2. For all θ ∈ Θ, the covariance function K θ is stationary. For all fixed t ∈ R d , K θ (t) is p + 1 times continuously differentiable with respect to θ. For all i 1 , ..., i p ∈ N so that i 1 + ... + i p ≤ p + 1, there exists A i1,...,ip < +∞ so that for all t ∈ R d , θ ∈ Θ,

∂ i1 ∂θ i1 1 ... ∂ ip ∂θ ip p K θ (t) ≤ A i1,...,ip 1 + |t| d+1 .
There exists a constant C inf > 0 so that, for any θ ∈ Θ, δ θ ≥ C inf . Furthermore, δ θ is p + 1 times continuously differentiable with respect to θ. For all i 1 , ..., i p ∈ N so that i 1 + ... + i p ≤ p + 1, there exists B i1,...,ip < +∞ so that for all θ ∈ Θ,

∂ i1 ∂θ i1 1 ... ∂ ip ∂θ ip p δ θ ≤ B i1,...,ip .
In Condition 3.2, we require a differentiability order of p + 1 for K θ and δ θ with respect to θ. In the related context of [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF], where a well-specified covariance model is studied, consistency of ML and CV can be proved with a differentiability order of 1 only. [One can check that the proofs of Propositions 3.1 and 3.4 in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] require only the first order partial derivatives of the Likelihood function.] The reason for this difference is that, as discussed after Theorem 3.4, an additional technical difficulty is present here, compared to [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF]. The specific approach we use requires the condition of differentiability order of p + 1 and we leave open the question of relaxing it. Note, anyway, that many parametric covariance models are infinitely differentiable with respect to the covariance parameters, especially the Matérn model. In Condition 3.2, assuming that the covariance function and its derivatives vanish with distance with a polynomial rate of order d + 1 is not restrictive. Indeed, many covariance functions vanish at least exponentially fast with distance.

Finally, the condition that the noise variance δ θ is lower bounded uniformly in θ is crucial for our proof methods. Since we address here noisy observations of Gaussian processes, this condition is reasonable so that the results of Section 4 can cover a large variety of practical situations, some of which are listed in Section 2. Note that even when the Gaussian process under consideration is observed exactly, it can be desirable to incorporate an instrumental positive term δ θ in the parametric model, for numerical reasons or for not interpolating exactly the observed values [START_REF] Andrianakis | The effect of the nugget on Gaussian process emulators of computer models[END_REF]. Thus, Condition 3.2 could also be considered for Gaussian processes that are observed without noise.

Maximum Likelihood

In this paper, the analysis of the ML estimator in the misspecified case is based on the Kullback-Leibler divergence of the distribution of y assumed under (K θ , δ θ ), for θ ∈ Θ, from the true distribution of y. More precisely, conditionally to X, y has a N (0, R 0 ) distribution and is assumed to have a N (0, R θ ) distribution. The Kullback-Leibler divergence of the latter distribution from the former is, after multiplication by 2/n,

D n,θ := 1 n log det R θ R -1 0 + Tr R 0 R -1 θ -1. (4)
The normalized Kullback-Leibler divergence in ( 4) is equal to 0 if and only if R θ = R 0 and is strictly positive otherwise. It is interpreted as an error criterion for using (K θ , δ θ ) instead of (K 0 , δ 0 ), when making inference on the Gaussian process Y .

Note that D n,θ is here appropriately scaled so that, if for a fixed θ (K θ , δ θ ) = (K 0 , δ 0 ), D n,θ should generally not vanish, nor diverge to infinity under increasing-domain asymptotics. This can be shown for instance in the framework of [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF], by using the methods employed there. It is also well-known that, in the case of a regular grid of observation points for d = 1, D n,θ converges to a finite limit as n → +∞ [START_REF] Azencott | Series of Irregular Observations: Forecasting and Model Building[END_REF]). This limit is twice the asymptotic Kullback information in [START_REF] Azencott | Series of Irregular Observations: Forecasting and Model Building[END_REF] and is positive if (K θ (t), δ θ ) differs from (K 0 (t), δ 0 ) for at least one point t in the regular grid of observation points. Similarly, in the spatial sampling framework of Condition 2.3, we observe in the Monte Carlo simulations of Section 4 that the order of magnitude of (4) does not change when n increases, for

(K θ , δ θ ) = (K 0 , δ 0 ).
The following theorem shows that the ML estimator asymptotically minimizes the normalized Kullback-Leibler divergence.

Theorem 3.3. Under Conditions 2.3, 3.1 and 3.2, we have, as n → ∞,

D n, θML = inf θ∈Θ D n,θ + o p (1),
where the o p (1) in the above display is a function of X and y only that goes to 0 in probability as n → ∞.

Theorem 3.3 is in line with the well-known fact that, in the iid setting, ML asymptotically minimizes the Kullback-Leibler divergence (than does not depend on sample size) from the true distribution, within a misspecified parametric model [START_REF] White | Maximum likelihood estimation of misspecified models[END_REF]. Theorem 3.3 conveys a similar message, with the normalized Kullback-Leibler divergence that depends on the spatial sampling. As discussed above, the infimum in Theorem 3.3 is typically lower bounded as n → ∞ in the misspecified case.

Note that Theorem 3.3 can be shown, in increasing-domain asymptotics, under other spatial samplings than that of Condition 2.3 (e.g. for the randomly perturbed regular grid of [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF]). Nevertheless, to the best of our knowledge, in the context of Condition 2.3, Theorem 3.3 is not a simple consequence of the existing literature, and an original proof is provided in Section A.2.

An important message of this paper is that the normalized Kullback-Leibler divergence is one error criterion among others and may not be given a particular priority. Indeed, most importantly, this criterion addresses the distribution of the Gaussian process only at the observation points. Hence, the minimality of D n, θML gives no information on the optimality of the inference of the values of Y at new points, obtained from (K θML , δ θML ). Especially, θML is typically not optimal for the prediction error at new points, as is illustrated in Section 4.

Cross Validation

Let ŷθ (t) = E θ|X (Y (t)|y). With the n × 1 vector r θ (t) so that (r θ (t)) j = K θ (t -X j ), we have ŷθ (t) = r t θ (t)R -1 θ y.
Then, define the family of random variables

E n,θ = 1 n [0,n 1/d ] d (ŷ θ (t) -Y (t)) 2 dt, (5)
where the integral is defined in the L 2 sense since K 0 is continuous. We call the criterion (5) the integrated square prediction error. It is natural to consider that the first objective of the CV estimator θCV is to yield a small E n, θCV . If the observation points X 1 , ..., X n are regularly spaced, then this objective might however not be fulfilled. Indeed, the principle of CV does not really have grounds in this case, since the LOO prediction errors are not representative of actual prediction errors for new points. This fact is only natural and has been noted in e.g. [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF] and [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF]. If however the observation points X 1 , ..., X n are not regularly spaced, then it is shown numerically in [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF] that the CV estimator θCV can yield a small E n, θCV and, especially, smaller than E n, θML . The following theorem, which is the main contribution of this paper, supports this conclusion under increasing-domain asymptotics.

Theorem 3.4. Under Conditions 2.3, 3.1 and 3.2, we have, as n → ∞,

E n, θCV = inf θ∈Θ E n,θ + o p (1),
where the o p (1) in the above display is a function of X, y and Y only that goes to 0 in probability as n → ∞.

Before discussing further the content of Theorem 3.4, we say a few words about its proof ( Section A.3). The main difference compared to [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] is that, because the observation points are independent, there is no minimum distance between two different observation points and clusters of closely spaced observation points may appear. As a consequence, the maximum eigenvalues of, say, the matrices R θ , are not upper bounded, even in probability. To the best of our knowledge, no increasing-domain asymptotic results for Gaussian processes are yet available in the literature that both address spatial sampling without minimal distance between two different observation points and require only assumptions on the covariance functions, as Conditions 3.1 and 3.2 do. Especially, the references [START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF]; Cressie andLahiri (1993, 1996) work under non-trivial assumptions on the covariance matrices involved, and show that these assumptions are fulfilled for examples of spatial samplings for which the minimal distance between two different observation points is bounded away from zero. More recently, also [START_REF] Zheng | On the asymptotics of maximum likelihood estimation for spatial linear models on a lattice[END_REF] assumes that the largest eigenvalues of the covariance matrices in the parametric model are bounded as n → ∞. In [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF], the working assumptions concern only the covariance model, but there is also a minimum distance between two different observation points. Note finally that in [START_REF] Lahiri | Central limit theorems for weighted sums of a spatial process under a class of stochastic and fixed designs[END_REF]; [START_REF] Lahiri | Asymptotic distributions of M-estimators in a spatial regression model under some fixed and stochastic spatial sampling designs[END_REF]; [START_REF] Lahiri | Resampling methods for spatial regression models under a class of stochastic designs[END_REF], random spatial sampling with independent observation points is addressed, under increasing and mixed increasing-domain asymptotics.

Nevertheless, these references do not address covariance parameter estimation and do not deal with M-estimation where the criterion function involves n × n random matrices.

As a consequence, the proof we propose for Theorem 3.4 is original and we do not address the asymptotic distribution of the ML and CV estimators. We leave this problem open to further research. Note nevertheless that, in the misspecified case addressed here, the fact that the ML and CV estimators minimize two different criteria and are thus typically asymptotically different is, in our opinion, at least as important as their asymptotic distributions.

An important element in the proof of Theorem 3.4 is that the variable t in the expression of the integrated square prediction error E n,θ in (5) plays the same role as a new point X n+1 , uniformly distributed on [0, n 1/d ] d and independent of (X 1 , ..., X n ). Hence, using the symmetry of X 1 , ..., X n+1 , for fixed θ, the mean value of E n,θ is equal to the mean value of a modification of the CV criterion CV θ in (2), where there are n + 1 observation points instead of n. Thus, one can indeed expect that the CV estimator minimizing CV θ also asymptotically minimizes E n,θ . [The challenging part for proving Theorem 3.4 is to control the deviations of the criteria E n,θ and CV θ from their mean values, uniformly in θ.] This discussion is exactly the paradigm of CV, that uses the LOO errors as empirical versions of the actual prediction errors. On the other hand, if the observation points constitute for instance a regular grid, then the variable t in E n,θ has close to nothing in common with them, so that Theorem 3.4 would generally not hold. This stresses that CV is generally not efficient for regular sampling of observation points, as discussed above.

In ( 5), we stress that E n,θ and the observation vector y are defined with respect to the same Gaussian process Y . Thus, Theorem 3.4 gives a guarantee for the estimator θCV relatively to the predictions it yields for the actual Gaussian process at hand. Theorems 3.3 and 3.4, when considered together with the results in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF], draw the following conclusion, on comparing ML and CV for covariance parameter estimation: In the well-specified case, these two estimators can be compared under the common criterion of estimation error and ML is preferable over CV. On the other hand, in the misspecified case, there is not a unique comparison criterion. ML is asymptotically optimal for the normalized Kullback-Leibler divergence (which does not address prediction of the Gaussian process Y at new points) while, under the random sampling addressed here, CV is more efficient than ML for the prediction mean square error. This last finding also supports some references listed in Section 1 addressing specific situations similar to the misspecified case.

Beyond a comparison of ML and CV, Theorems 3.3 and 3.4 show that both estimators can not do significantly worse asymptotically than any other potential estimators for their respective error criteria. This last fact holds for ML under fairly general samplings of observation points, as we have discussed, but necessitates here a purely random sampling for CV.

Finally, the fact that ML and CV typically optimize different criteria in the misspecified case can serve as a practical guideline. That is, one can compute the estimated covariance parameters with both methods and compare the two estimates and the corresponding log-likelihood and LOO mean square error values. If the differences between ML and CV are large, then it could be a warning that the covariance model at hand can be inappropriate.

Monte Carlo simulation

We illustrate Theorems 3.3 and 3.4 in a Monte Carlo simulation. We consider the Matérn covariance model in dimension d = 1. A covariance function on R is Matérn (σ 2 , , ν) when it is written

K σ 2 , ,ν (t) = 1 Γ (ν) 2 ν-1 2 √ ν |t| ν K ν 2 √ ν |t| ,
with Γ the Gamma function and K ν the modified Bessel function of second order. The parameters σ 2 , and ν are respectively the variance, correlation length and smoothness parameters. We refer to e.g [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] for a presentation of the Matérn covariance function.

In the simulation, the true covariance function of Y is Matérn (σ 2 0 , 0 , ν 0 ) with σ 2 0 = 1, 0 = 3 and ν 0 = 10. This choice of ν 0 corresponds to a smooth Gaussian process and enables, as we see below, to illustrate Theorems 3.3 and 3.4 in a more striking manner. The true noise variance is δ 0 = 0.25 2 .

For the covariance and noise variance model, it is considered that the smoothness parameter ν 0 is known, that the noise variance δ 1 is fixed, and that the parameter θ = (σ 2 , ) is estimated by ML or CV. For both ML and CV, the optimization is restricted to the domain Θ = [0.1 2 , 10 2 ] × [0.2, 10]. [We experience that the conclusions of the Monte Carlo simulation are the same if a larger optimization domain is considered.] The well-specified case corresponds to δ 1 = δ 0 and the misspecified case corresponds to δ 1 = 0.1 2 = δ 0 . This covariance model is representative of practical applications. Indeed, first it is common practice to fix the value of the smoothness parameter in the Matérn model, as is discussed in Section 2. Second, when using Gaussian process models on experimental or natural data, it can often occur that field experts provide an a priori value for the noise variance (see e.g. [START_REF] Bachoc | Calibration and improved prediction of computer models by universal Kriging[END_REF]). The misspecified case we address corresponds to an underestimation of the noise variance, possibly because some sources of measurement errors have been neglected.

The Monte Carlo simulation is carried out as follows. For n = 100 and N = 1000 or n = 500 and N = 200 we repeat N data generations, estimations and quality criterion computations and average the results. More specifically, we simulate N independent realizations of the observation points, of the observation vector and of the Gaussian process on [0, n], under the true covariance function and noise variance. For each of these N realizations, we compute the ML and CV estimates under the well-specified and misspecified models. For each of these estimates of the form θ = (σ 2 , ˆ ), we compute the corresponding criteria D n,σ 2 , ˆ and E n,σ 2 , ˆ .

In Figure 1 we report, for n = 100, for the well-specified and misspecified cases and for ML and CV, the histograms of the estimates ˆ , and of the values of the error criteria D n,σ 2 , ˆ and E n,σ 2 , ˆ . In the well-specified case, the conclusions are in agreement with the main message of previous literature: Both estimators estimate the true 0 = 3 with reasonable accuracy and have error criteria that are relatively small. By also considering Table 1, where the averages and standard deviations corresponding to Figures 1 and2 are reported, we observe that ML performs better than CV in all aspects. The estimation error for and the normalized Kullback-Leibler divergence are significantly smaller for ML, while the integrated square prediction error is similar under ML and CV estimation, but nonetheless smaller for ML.

The conclusions are however radically different in the misspecified case, as is implied by Theorems 3.3 and 3.4. First, the ML estimates of are significantly smaller than in the well-specified case, and can even be equal to the lower-bound 0.2. The ML estimates of σ 2 are not reported in Figure 1 to save space and are close to 1, so that, approximately, the variance of the observations, as estimated by ML, is close to the true variance of the observations. The reason for these small estimates of by ML is the underestimation of the noise variance δ 0 , coupled with the large smoothness parameter ν 0 . Indeed, there exist pairs of closely spaced observation points for which the corresponding differences of observed values are large compared to δ 1 , so that for values of that are larger than those computer by ML, the criterion (1) blows up, for all values of σ 2 . [Using a value of σ 2 smaller or approximately equal to 1 does not counterbalance the damaging impact on (1) of these pairs of closely spaced observation points with large observed value differences. Increasing σ 2 over 1 is also not optimal for (1), since on a large scale, the observations do have variances close to 1.] This phenomenon for ML is all the more important when the smoothness parameter ν 0 is large, which is why we choose here the value ν 0 = 10 to illustrate it. To summarize, ML gives an important weight to pairs of closely spaced observation points with large observation differences and consequently estimates small correlation lengths to explain, so to speak, these observation differences.

On the contrary for CV, if we consider only the predictions ŷσ 2 , (t) at new points t and the LOO predictions ŷi,σ 2 , , with ( , σ 2 ) ∈ Θ, then the situation is virtually the same as if the model was well-specified. Indeed, the covariance matrices and vectors obtained from σ 2 , and δ 0 are equal to δ 0 /δ 1 time those obtained from σ 2 δ 1 /δ 0 , and δ 1 , so that the corresponding predictions are identical. Hence, in Figure 1, the empirical distribution of ˆ CV is approximately the same between the well-specified and misspecified cases. In the misspecified case, we find that the empirical distribution of σ2

CV (not reported in Figure 1 to save space) is δ 1 /δ 0 time that of the well-specified case. Of course, although the CV predictions are not damaged by the misspecified δ 1 , the CV estimations of other characteristics of the conditional distribution of Y given the observed data are damaged.

[For example the confidence intervals for Y (t) obtained from the CV estimates are significantly too small.]

The histograms of E n,σ 2 , ˆ and D n,σ 2 , ˆ for ML and CV in Figure 1 confirm the discussion on the estimated parameters. For ML which estimates small correlation lengths, the error criteria E n,σ 2 M L , ˆ M L are significantly larger than in the well-specified case (by an approximate factor of 3 on average as seen in Table 1). The error criteria D n,σ 2 M L , ˆ M L also increase and become larger than these of both ML and CV in the well-specified case. For CV, the error criteria E n,σ 2 CV , ˆ CV are, as discussed, as small as in the well-specified case and approximately 3 times smaller on average than for ML, illustrating Theorem 3.4. However, the error criteria D n,σ 2 CV , ˆ CV are 3 times larger for CV than for ML, in the misspecified case, illustrating Theorem 3.3.

Finally, in Figure 2, the settings are the same as for Figure 1 but for n = 500. The conclusions on the comparison between ML and CV are the same as for n = 100. Then, the estimates of under ML and CV have less variance than for n = 100, and their histograms are approximately unimodal and symmetric. Finally, for ML and CV in the misspecified case, E n,σ 2 , ˆ and D n,σ 2 , ˆ keep the same averages between n = 100 and n = 500. In the well-specified case, E n,σ 2 , ˆ also keeps the same average, while D n,σ 2 , ˆ becomes very small. This is because D n,σ 2 0 , 0 = 0 in the well-specified case, while E n,σ 2 0 , 0 is non-zero and should not vanish to 0 as n → ∞, since the density of observation points in the prediction domain is constant with n.

Discussion

We prove that for independent and uniformly distributed observation points, CV asymptotically yields the smallest integrated square prediction error among all possible parameters. On the contrary ML is asymptotically optimal for the Kullback-Leibler divergence which does not provide information for prediction at new points. In the case addressed in the Monte Carlo simulation, it is clear that the covariance parameters estimated by ML and CV, and their corresponding performances for the integrated square prediction error and the Kullback-Leibler divergence are significantly different. Thus, the contrast between the well-specified case, where both estimators converge to the true covariance parameter, is striking.

For regular sampling, the principle of CV does not really have grounds while ML would typically still asymptotically minimize the Kullback-Leibler divergence. Note however that regular sampling can still damage the ML estimation and make model misspecification more difficult to detect (see [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of gaussian processes[END_REF] and the references discussed in its introduction and conclusion).

Finally, the results of the Monte Carlo simulation make it conceivable that, for independent and uniform observation points, the ML and CV estimators converge to optimal parameters, for respectively the Kullback-Leibler divergence and the integrated square prediction error, and are asymptotically normal. [These optimal parameters would be equal to the true ones in the well-specified case.] Considering asymptotic normality might require new techniques to account for the absence of minimum distance between different observation points. We leave this problem open to future research.

A Proofs

A.1 Notation

In all the appendix, we consider that Conditions 2.3, 3.1 and 3. of the Gaussian process Y on [0, n] with Matérn (σ 2 0 = 1, 0 = 3, ν 0 = 10) covariance function, and of the corresponding observation vector with noise variance δ 0 = 0.25 2 . For each simulation, ν 0 is known, the noise variance is fixed to δ 1 = δ 0 (well-specified case) or δ 1 = 0.1 2 = δ 0 (misspecified case), σ 2 and are estimated by ML and CV and the corresponding error criteria D n,σ 2 , ˆ (normalized Kullback-Leibler divergence) and E n,σ 2 , ˆ (integrated square prediction error) are computed. The histograms of the N estimates of and of the N corresponding values of the error criteria are reported for ML and CV and in the well-specified and misspecified cases. In the well-specified case, the estimates are on average reasonably close to the true values, the error criteria are reasonably small and ML performs better than CV in all aspects. In the misspecified case, the ML and CV estimates of the correlation lengths are significantly different, ML performs better than CV for D n,σ 2 , ˆ and CV performs better than ML for E n,σ 2 , ˆ . Between Figure 1 and Figure 2, all error criteria keep the same order of magnitude except D n,σ 2 , ˆ in the well-specified case which decreases and becomes very small for ML and CV (see also Table 1). This is because in the well-specified case D n,σ 2 0 , 0 = 0. For a sequence of real random variables z n , we write z n → p 0 and z n = o p (1) when z n converges to zero in probability. For a random variable A and a deterministic function f (A), we may write E A (f (A)) for E(f (A)). For two random variables A and B and a deterministic function f (A, B) we may write E A|B (f (A, B)) for E(f (A, B)|B).

|A| 2 = 1 m m i,j=1 A 2 i,
For a finite set E, we write |E| for its cardinality. For a continuous set E ⊂ R d , we write |E| for its Lebesgue measure. For two sets A, B in R d , we write d(A, B) = inf a∈A,b∈B |a -b|.

We write C sup a generic non-negative finite constant (not depending on n, X, Y , and θ). The actual value of C sup is of no interest and can change in the same sequence of equations. For instance, instead of writing, say, a ≤ 2b ≤ 4c, we shall write a ≤ C sup b ≤ C sup c. Similarly, we write C inf a generic strictly positive constant (not depending on n, X, Y , and θ).

A.2 Proofs for Maximum Likelihood

Lemma A.1. For i = 1, ..., p,

E sup θ∈Θ ∂ ∂θ i L θ is bounded w.r.t n.
Proof of Lemma A.1. [START_REF] Gray | Toeplitz and circulant matrices: A review[END_REF] and R -1

E sup θ∈Θ ∂ ∂θ i L θ = E sup θ∈Θ 1 n Tr R -1 θ ∂ ∂θ i R θ - 1 n y t R -1 θ ∂ ∂θ i R θ R -1 θ y (Cauchy-Schwarz:) ≤ E   sup θ∈Θ R -1 θ 2 ∂ ∂θ i R θ 2   (6) + 1 n E sup θ∈Θ R -1 θ y 2 1 n E sup θ∈Θ ∂ ∂θ i R θ R -1 θ y 2 (7) Now, R -1 θ 2 ≤ R -1 θ 2 because of (2.19) in
θ 2 is bounded uniformly in θ because of Lemma A.18. Also, E sup θ∈Θ |[∂/∂θ i ]R θ |
2 is bounded because of Condition 3.2 and of a simple case of Lemma A.17. So the right-hand side of ( 6) is bounded because of Jensen inequality. It remains to show that the term ( 7) is bounded. To show this, note first that

1 n E sup θ∈Θ R -1 θ y 2 ≤ 1 n E sup θ∈Θ ||y|| 2 R -1 θ 2 (Lemma A.18:) C sup n E ||y|| 2 = C sup (K 0 (0) + δ 0 ), is bounded. Thus, it remains to show that 1 n E sup θ∈Θ ([∂/∂θ i ]R θ ) R -1 θ y 2 is bounded.
For this, we have

1 n E sup θ∈Θ ∂ ∂θ i R θ R -1 θ y 2 = 1 n E sup θ∈Θ y t R -1 θ ∂ ∂θ i R θ 2 R -1 θ y ≤ C sup i1+...+ip≤p Θ 1 n E ∂ i1 ∂θ i1 1 ... ∂ ip ∂θ ip p y t R -1 θ ∂ ∂θ i R θ 2 R -1 θ y dθ (Lemma A.15).
Thus, it suffices to show that, for fixed i 1 , ..., i p ∈ N so that i 1 + ... + i p ≤ p,

1 n sup θ∈Θ E ∂ i1 ∂θ i1 1 ... ∂ ip ∂θ ip p y t R -1 θ ∂ ∂θ i R θ 2 R -1 θ y
is bounded. The above display is smaller than a fixed sum of terms of the form (1/n) sup θ∈Θ E (|y t M θ y|), where the number of terms is independent of n and M θ is of the form

N 1,θ M 1,θ ...M k,θ N k+1,θ with N i,θ = I n or N i,θ = R -1 θ and with M i,θ of the form [∂ c1 /∂θ c1 1 ]...[∂ cp /∂θ cp p ]R θ with c 1 , .
.., c p ∈ N and c 1 + ... + c p ≤ p + 1. Hence, it is enough to show that any term of the form sup θ∈Θ (1/n)E (|y t M θ y|) above is bounded. We have

sup θ∈Θ 1 n E |y t M θ y| ≤ sup θ∈Θ 1 n E y t M θ y -E y t M θ y X + sup θ∈Θ 1 n E E y t M θ y X ≤ sup θ∈Θ E var 1 n y t M θ y X + sup θ∈Θ 1 n E E y t M θ y X (Jensen inequality) = sup θ∈Θ E 1 2n 2 Tr [R 0 {M θ + M t θ } R 0 {M θ + M t θ }] + sup θ∈Θ 1 n E (|Tr [R 0 M θ ]|) ≤ sup θ∈Θ 1 2n E |R 0 {M θ + M t θ }| 2 + sup θ∈Θ E (|R 0 | 2 ) E (|M θ | 2 ) (Cauchy-Schwarz). ≤ sup θ∈Θ 1 n E |R 0 M θ | 2 + |R 0 M t θ | 2 + sup θ∈Θ E (|R 0 | 2 ) E (|M θ | 2 )
In the display above, the first term goes to 0 because of Conditions 3.1 and 3.2 and Lemmas A.17 and A.18. The second term is bounded because of Lemmas A.17, A.18 and A.21. This completes the proof.

Corollary A.2. For any i = 1, ..., p,

E sup θ∈Θ ∂ ∂θ i E ( L θ | X) and sup θ∈Θ ∂ ∂θ i D n,θ are bounded w.r.t n. Proof of corollary A.2. Note first that [∂/∂θ i ]E ( L θ | X) = [∂/∂θ i ]D n,θ .
The corollary is then a consequence of the fact that, for fixed n, we have

(∂/∂θ i )E ( L θ | X) = E ( (∂/∂θ i )L θ | X) and of sup θ |E(.)| ≤ E(sup θ |.|). Lemma A.3. Consider a fixed θ ∈ Θ. Then E (|L θ -E(L θ |X)|) → n→∞ 0.
Proof of Lemma A.3. We have, applying Jensen inequality twice

E (|L θ -E(L θ |X)|) ≤ E var(L θ |X) ≤ E (var(L θ |X)) = E 2 n 2 Tr R 0 R -1 θ R 0 R -1 θ .
The eigenvalues of R -1 θ are smaller than a finite constant C sup for any n, X, θ from Lemma A.18. Thus, by applying Cauchy schwarz inequality and Lemmas A.17 and A.21,

E (|L θ -E(L θ |X)|) ≤ 2 n E |R 0 R -1 θ | 2 ≤ C sup √ n .
Proof of Theorem 3.3. We have

sup θ∈Θ |L θ -log(det(R 0 )) -1 -D n,θ | ≤ sup θ∈Θ |L θ -E(L θ |X)| + sup θ∈Θ |E(L θ |X) -log(det(R 0 )) -1 -D n,θ | = sup θ∈Θ |L θ -E(L θ |X)| .
The term in the above display goes to 0 in probability. Indeed, for fixed θ, the function of θ goes to 0 in probability because of Lemma A.3. The convergence of the supremum over θ to 0 is then a consequence of the fact that Θ is compact and of Lemma A.1 and corollary A.2. Finally, since θML minimizes L θ and so also L θ -log(det(R 0 )) -1, we conclude with, for any θ ∈ Θ,

D n, θML -D n,θ ≤ L θML -L θ + 2 sup θ∈Θ |L θ -log(det(R 0 )) -1 -D n,θ | ≤ 2 sup θ∈Θ |L θ -log(det(R 0 )) -1 -D n,θ | .
Hence sup θ∈Θ D n, θML -D n,θ = o p (1).

A.3 Proofs for Cross Validation

Lemma A.4. For i = 1, ..., p,

E sup θ∈Θ ∂ ∂θ i CV θ is bounded w.r.t n.
Proof of Lemma A.4. We have

E sup θ∈Θ ∂ ∂θ i CV θ ≤ E 1 n n k=1 sup θ∈Θ ∂ ∂θ i (y k -ŷk,θ ) 2 (symmetry of X 1 , ..., X n :) = E sup θ∈Θ ∂ ∂θ i (y 1 -ŷ1,θ ) 2 (Lemma A.15:) ≤ C sup i1+...+ip≤p Θ E ∂ i1 ∂θ i1 1 ... ∂ ip ∂θ ip p ∂ ∂θ i (y 1 -ŷ1,θ ) 2 dθ. Let us consider a specific i 1 , ..., i p . Then [∂ i1 /∂θ i1 1 ]...[∂ ip /∂θ ip p ][∂/∂θ i ](y 1 -ŷ1,θ
) 2 is a weighted sum (weights and number of terms depending only on i 1 , ..., i p ), so that the terms are of the two following forms:

(y 1 -ŷ1,θ ) ∂ k1 ∂θ k1 1 ... ∂ kp ∂θ kp p ∂ ∂θ i ŷ1,θ or ∂ k1 ∂θ k1 1 ... ∂ kp ∂θ kp p ∂ ∂θ i ŷ1,θ ∂ l1 ∂θ l1 1 ... ∂ lp ∂θ lp p ∂ ∂θ i ŷ1,θ .
Thus, we just have to show that the mean values of the absolute values of the terms of the form above (for k 1 + ... + k p ≤ p and l 1 + ... + l p ≤ p) are bounded uniformly in θ ∈ Θ. By using Cauchy-Schwarz inequality, these means of absolute values are smaller than either

E ((y 1 -ŷ1,θ ) 2 ) E   ∂ k1 ∂θ k1 1 ... ∂ kp ∂θ kp p ∂ ∂θ i ŷ1,θ 2   or E   ∂ k1 ∂θ k1 1 ... ∂ kp ∂θ kp p ∂ ∂θ i ŷ1,θ 2   E   ∂ l1 ∂θ l1 1 ... ∂ lp ∂θ lp p ∂ ∂θ i ŷ1,θ 2   . Now, E (y 1 -ŷ1,θ ) 2 ≤ 2E(y 2 1 ) + 2E(ŷ 2 1,θ ).
The term E(y 2 1 ) is bounded uniformly in θ. Thus, finally, it remains to show that for any a 1 + ...

+ a p ≤ p + 1, sup θ∈Θ E [∂ a1 /∂θ a1 1 ]...[∂ ap /∂θ ap p ]ŷ 1,θ 2 is bounded. For that, we have ŷ1,θ = r t 1,θ R -1 1,θ y -1 . Thus, [∂ a1 /∂θ a1 1 ]...[∂ ap /∂θ ap p ]ŷ 1,θ is a fixed sum of weighted terms of the form w t θ M θ y -1 , where w θ is of the form [∂ b1 /∂θ b1 1 ]...[∂ bp /∂θ bp p ]r 1,θ (b 1 + ... + b p ≤ p + 1) and M θ is of the form R -1 1,θ M 1,θ ...R -1 1,θ M k,θ R -1 1,θ . Finally, k is smaller than a finite constant C sup (function of p) and M i,θ is of the form [∂ c1 /∂θ c1 1 ]...[∂ cp /∂θ cp p ]R 1,θ , with c 1 + ... + c p ≤ p + 1.
Thus, it is sufficient to show that a generic sup θ∈Θ E (w t θ M θ y -1 ) 2 , as previously defined, is bounded.

Then,

sup θ∈Θ E w t θ M θ y -1 2 = sup θ∈Θ E X E y|X y t -1 M t θ w θ w t θ M θ y -1 = sup θ∈Θ E X Tr R 1,0 M t θ w θ w t θ M θ ≤ sup θ∈Θ E X   n i,j=2 M θ R 1,0 M t θ i,j w θ w t θ i,j   (8) = sup θ∈Θ   n i,j=2 E X2,...,Xn M θ R 1,0 M t θ i,j E X1|X2,...,Xn w θ w t θ i,j   .
Now, because of Conditions 2.3 and 3.2,

E X1|X2,...,Xn w θ w t θ i,j ≤ C sup n [0,n 1/d ] d 1 1 + |X i -x 1 | d+1 1 1 + |X j -x 1 | d+1 dx 1 (Lemma A.16:) ≤ 1 n C sup 1 + |X i -X j | d+1 . So, sup θ∈Θ E w t θ M θ y -1 2 ≤ C sup 1 n sup θ∈Θ   n i,j=2 E X2,...,Xn M θ R 1,0 M t θ i,j 1 1 + |X i -X j | d+1   (Cauchy-Schwarz:) ≤ sup θ∈Θ    E |M θ R 1,0 M t θ | 2 E    1 n n i,j=2 1 1 + |X i -X j | d+1 2       .
The supremum over θ of the second term above is bounded because of Lemma A.17. The supremum over θ of the first term above is bounded because of Lemmas A.17, A.18 and A.21.

Corollary A.5. For any i = 1, ..., p,

E sup θ∈Θ ∂ ∂θ i E ( CV θ | X) and sup θ∈Θ ∂ ∂θ i E (CV θ ) are bounded w.r.t n.
Proof of corollary A.5. The corollary is a consequence of Lemma A.4, sup θ |E(.)| ≤ E(sup θ |.|) and of the fact that, for fixed n, we have

(∂/∂θ i )E ( CV θ | X) = E ( (∂/∂θ i )CV θ | X) and (∂/∂θ i )E(CV θ ) = E((∂/∂θ i )CV θ ).
Lemma A.6. For any fixed θ ∈ Θ we have

E (|CV θ -E(CV θ |X)|) → n→∞ 0.
Proof of Lemma A.6. From (3), we have

CV θ = y t M θ y, with M θ = (1/n)R -1 θ diag(R -1 θ ) -2 R -1 θ .
Because of Lemma A.18, the eigenvalues of M θ are bounded uniformly in n, X, θ by a finite constant C sup . Thus, the proof of the lemma is exactly the same as that of Lemma A.3, with R -1 θ replaced by M θ .

Definition A.7. Consider a fixed θ ∈ Θ. Consider two functions of n: n 2 (n) ∈ N * and ∆(n) ≥ 0, that we write n 2 and ∆ for simplicity, so that, for any n ∈ N * , n 2 can be written n 2 = N d 2 , with N 2 ∈ N * , and so that n = n 2 ∆. Let, for i = 1, ..., N 2 -1, =C(t2) . Define Rθ , Ri,θ , ri,θ , ỹi,θ , CV θ similarly to R θ , R i,θ , r i,θ , ŷi,θ , CV θ but with K θ replaced by Kθ . Furthermore, let us write the n 2 aforementioned sets of the form d j=1 c ij , for i 1 , ..., i d ∈ {1, ..., N 2 }, as the sets C 1 , ..., C n2 . [The specific one-to-one correspondence we use between {1, ..., N 2 } d and {1, ..., n 2 } is of no interest. Note that this one-to-one correspondence depends on n. The sets C 1 , ..., C n2 also depend of n, but we drop this dependence in the notation for simplicity.]

c i = [((i -1)/N 2 )n 1/d , (i/N 2 )n 1/d ). Let c N2 = [((N 2 -1)/N 2 )n 1/d , n 1/d ]. Let, for x ∈ [0, n 1/d ], i(x) be the unique i ∈ {1, ..., N 2 } so that x ∈ c i . Let, for t = (t 1 , ..., t d ) t ∈ [0, n 1/d ] d , C(t) = d j=1 c i(tj ) . Define the non-stationary covariance function Kθ (t 1 , t 2 ) = K θ (t 1 , t 2 )1 C(t1)
Let N i be the random number of observation points in C i and let X i be the random N i -tuple obtained from X by keeping only the observation points that are in C i and by preserving the order of the indices in X. Let y i be the column vector of size N i , composed by the components y j of y for which X j is in C i (preserving the order of indexes). Let Ri,θ and Ri,0 be the covariance matrices, under (K θ , δ θ ) and (K 0 , δ 0 ), of y i , given X.

Finally, for 1 ≤ i, j ≤ n 2 , let v i and w j be two N i × 1 and N j × 1 vectors and M ij be a N i × N j matrix. Then we use the convention that, when

N i = 0, |M ij | = ||M ij || = 0, ||v i || = |v i | = 0 and v t i M ij w j = 0. Furthermore, if i = j and M ii is invertible when N i ≥ 1, we use the convention that v t i (M ii ) -1 w i = 0 when N i = 0.
[These conventions enable to write equalities or inequalities involving matrices and vectors of size N i , N j or N i × N j , that hold regardless of whether N i or N j are zero or not. As can be checked along the proofs involving Definition A.7, these relations boil down to trivial relations (e.g. 0 = 0) when N i = 0 or N j = 0. This way of proceeding considerably simplifies the exposition in these proofs.] Lemma A.8. Consider a fixed θ ∈ Θ. In the context of Definition A.7, if n 2 = o(n),

E CV θ -CV θ → n→∞ 0.
Proof of Lemma A.8. Assume that n 2 = o(n), or equivalently that ∆ → n→∞ ∞. We have

E CV θ -CV θ ≤ 1 n n i=1 E (y i -ŷi,θ ) 2 -(y i -ỹi,θ ) 2 (symmetry:) = E (y 1 -ŷ1,θ ) 2 -(y 1 -ỹ1,θ ) 2 = E (y 1 -r t 1,θ R -1 1,θ y -1 ) 2 -(y 1 -rt 1,θ R-1 1,θ y -1 ) 2 = E rt 1,θ R-1 1,θ y -1 -r t 1,θ R -1 1,θ y -1 2y 1 -r t 1,θ R -1 1,θ y -1 -rt 1,θ R-1 1,θ y -1 (Cauchy-Schwarz:) ≤ E rt 1,θ R-1 1,θ y -1 -r t 1,θ R -1 1,θ y -1 2 E 2y 1 -r t 1,θ R -1 1,θ y -1 -rt 1,θ R-1 1,θ y -1 2
Now, the second square root in the above display is bounded, because of (a + b + c) 2 ≤ 3 a 2 + b 2 + c 2 and of arguments similar to but simpler than those given in the proof of Lemma A.4. Thus it only remains to show that

E rt 1,θ R-1 1,θ y -1 -r t 1,θ R -1 1,θ y -1 2 → n→∞ 0.
For this,

E rt 1,θ R-1 1,θ y -1 -r t 1,θ R -1 1,θ y -1 2 ≤ 2E rt 1,θ ( R-1 1,θ -R -1 1,θ )y -1 2 (9) +2E (r 1,θ -r 1,θ ) t R -1 1,θ y -1 2 .
We show separately that both terms in the right-hand side of (9) converge to 0. For the first term,

E rt 1,θ ( R-1 1,θ -R -1 1,θ )y -1 2 = E Tr R 1,0 ( R-1 1,θ -R -1 1,θ )r 1,θ rt 1,θ ( R-1 1,θ -R -1 1,θ ) ≤ n i,j=2 E |( R-1 1,θ -R -1 1,θ )R 1,0 ( R-1 1,θ -R -1 1,θ )| i,j |r 1,θ rt 1,θ | i,j .
Hence, by the same arguments as after (8) in the proof of Lemma A.4, we obtain

E rt 1,θ ( R-1 1,θ -R -1 1,θ )y -1 2 2 ≤ C sup E |( R-1 1,θ -R -1 1,θ )R 1,0 ( R-1 1,θ -R -1 1,θ )| 2 ≤ C sup E || R-1 1,θ || + ||R -1 1,θ || |R 1,0 ( R-1 1,θ -R -1 1,θ )| 2 (Lemmas A.18 and A.19:) ≤ C sup n E Tr ( R-1 1,θ -R -1 1,θ )R 2 1,0 ( R-1 1,θ -R -1 1,θ ) (Cauchy-Schwarz:) ≤ C sup E |( R-1 1,θ -R -1 1,θ ) 2 | 2 E |R 2 1,0 | 2 . From Lemma A.17, E |R 2 1,0 | 2 is bounded, so it remains to show that E |( R-1 1,θ -R -1 1,θ ) 2 | 2 converges to 0. For this, E |( R-1 1,θ -R -1 1,θ ) 2 | 2 = E |( R-1 1,θ (R 1,θ -R1,θ )R -1 1,θ ) 2 | 2 (Lemma A.18:) ≤ C sup E |(R 1,θ -R1,θ )R -1 1,θ R-1 1,θ (R 1,θ -R1,θ )| 2 = C sup 1 n E Tr (R 1,θ -R1,θ ) 2 R -1 1,θ R-1 1,θ (R 1,θ -R1,θ ) 2 R-1 1,θ R -1 1,θ (Cauchy-Schwarz:) ≤ C sup E |(R 1,θ -R1,θ ) 2 R -1 1,θ R-1 1,θ | 2 E |(R 1,θ -R1,θ ) 2 R-1 1,θ R -1 1,θ | 2 .
Hence, with Lemmas A.18, A.19 and A.22, we conclude that the first term of the right hand side of (9) goes to 0. Let us now show that the second term of the right hand side of (9) goes to 0. We have,

E (r 1,θ -r 1,θ ) t R -1 1,θ y -1 2 = E Tr R -1 1,θ R 1,0 R -1 1,θ (r 1,θ -r 1,θ )(r 1,θ -r 1,θ ) t ≤ n i,j=2 E [R -1 1,θ R 1,0 R -1 1,θ ] i,j 1 n [0,n 1/d ] d 1 1 + |X i -x 1 | d+1 1 1 + |X j -x 1 | d+1 1 C(Xi) =C(x1) 1 C(Xj ) =C(x1) dx 1 ,
where the last line is obtained similarly to after (8) in the proof of Lemma A.4. Thus we have, with the notation and result of Lemma A.23,

E (r 1,θ -r 1,θ ) t R -1 1,θ y -1 2 ≤ C sup 1 n n i,j=2 E [R -1 1,θ R 1,0 R -1 1,θ ] i,j 1 1 + |X i -X j | d+1 f (D ∆ (X i , X j )) (Cauchy-Schwarz:) ≤ E |R -1 1,θ R 1,0 R -1 1,θ | 2 1 n n i,j=2 E 1 1 + |X i -X j | d+1 2 f 2 (D ∆ (X i , X j )) .
From Lemmas A.17, A.18 and A.21, the first √ . in the above display is bounded. Thus it remains to show that the second √ . goes to 0. For this, noting that f 2 (t) ≤ C sup f (t) and distinguishing the case i = j from the case i = j,

1 n n i,j=2 E 1 1 + |X i -X j | d+1 2 f (D ∆ (X 1 , X j )) ≤ C sup n [0,n 1/d ] d f (D ∆ (x))dx + C sup n [0,n 1/d ] d dx 1 [0,n 1/d ] d dx 2 1 1 + |x 1 -x 2 | d+1 f (D ∆ (x 1 , x 2 )) = C sup n [0,n 1/d ] d f (D ∆ (x))dx + o(1) (Lemma A.24). (10)
Now, for any > 0, there is a finite T so that f (T ) ≤ , and by defining

E n = {x ∈ [0, n 1/d ] d ; D ∆ (x) ≤ T }, we have |E n | = o(n),
as can be seen easily, and

1 n [0,n 1/d ] d f (D ∆ (x))dx ≤ f (0) |E n | n + .
This finally shows that the second term of the right hand side of (9) goes to 0 which finishes the proof.

Lemma A.9. For any fixed θ ∈ Θ,

E (|E (CV θ ) -E [CV θ |X]|)
goes to 0 as n → ∞.

Proof of Lemma A.9. Fix θ ∈ Θ. Because of Lemma A.8 and of |E(.)| ≤ E(|.|), it is sufficient to show that there exists a sequence ∆ → +∞ so that the lemma holds with CV θ replaced by CV θ . Then, because of (E(.))

2 ≤ E (.) 2 , it is sufficient to show var E CV θ X → n→∞ 0.

Let C 1 , ..., C n2 be as in Definition A.7. Define, for k = 1, ..., n 2 ,

f k (X) = 1 ∆ Xi∈C k E y i -ỹi,θ 2 X .
[Note that, following the discussion in Definition A.7, we have

f k (X) = 0 if N k = 0 and f k (X) = K 0 (0) + δ 0 if N k = 1.] Then E CV θ |X = (1/n 2 ) n2 k=1 f k (X).
Let Rk,θ and Rk,0 be as in Definition A.7. Because of the definition of K and by (3), we have

f k (X) = 1 ∆ Tr Rk,0 R-1 k,θ diag( R-1 k,θ ) -2 R-1 k,θ .
The functions f k (X) satisfy the conditions of Lemma A.25. Furthermore, by using the notation N k of Lemma A.25, we have

E f 2 k (X) N k = N = 1 ∆ 2 E Tr Rk,0 R-1 k,θ diag( R-1 k,θ ) -2 R-1 k,θ 2 N k = N (Cauchy-Schwarz:) ≤ 1 ∆ 2 E N 2 | Rk,0 | 2 | R-1 k,θ diag( R-1 k,θ ) -2 R-1 k,θ | 2 N k = N (Lemma A.20:) ≤ C sup N 2 ∆ 2 E | Rk,0 | 2 N k = N (Condition 3.1 and Lemma A.25:) ≤ C sup N 2 ∆ 2 1 + N ∆ 2 [0,∆ 1/d ] d [0,∆ 1/d ] d 1 1 + |x 1 -x 2 | d+1 dx 1 dx 2 ≤ C sup N 2 ∆ 2 + N 3 ∆ 3 ≤ C sup 1 + N 4 ∆ 4 .
Thus, because of Lemma A.26, there exists a sequence ∆ → n→∞ ∞ so that var E CV θ X → n→∞ 0, which completes the proof.

Lemma A.10. Let, with the notation of Definition A.7, Ẽn,θ be defined as E n,θ , with K θ replaced by Kθ . Fix

θ ∈ Θ. Then, if n 2 = o(n), E E n,θ -Ẽn,θ → n→∞ 0.
Proof of Lemma A.10. We have, by letting ỹθ (t) be as ŷθ (t), with K θ replaced by Kθ .

E E n,θ -Ẽn,θ = E 1 n [0,n 1/d ] d [Y (t) -ŷθ (t)] 2 dt - 1 n [0,n 1/d ] d Y (t) -ỹθ (t) 2 dt ≤ E 1 n [0,n 1/d ] d [Y (t) -ŷθ (t)] 2 -Y (t) -ỹθ (t) 2 dt
The variable t in the integral above is formally equivalent to a new observation point X n+1 , so that X 1 , ..., X n+1 are independent and uniformly distributed on [0, n 1/d ] d . Thus,

E E n,θ -Ẽn,θ ≤ E (Y (X n+1 ) -ŷθ (X n+1 )) 2 -Y (X n+1 ) -ỹθ (X n+1 ) 2 .
The rest of the proof is carried out as in Lemma A.8, the only difference being that there are n + 1 observation points instead of n.

Lemma A.11. For any fixed θ ∈ Θ we have

E (|E n,θ -E(E n,θ |X)|) → n→∞ 0.
Proof of Lemma A.11. Because of Lemma A.10 and using |E(.)| ≤ E(|.|) and E 2 (.) ≤ E((.) 2 ), it is sufficient to show that there exists a sequence ∆ → n→∞ so that E var Ẽn,θ X goes to 0 as n → ∞. Now, using the notation N i of Lemma A.25 and Conditions 3.1 and 3.2,

T 2 ≤ C sup 1 n 2 n2 i=1 1 n 2 n2 j=1 1 ∆ 2 Ci dt i Cj dt j E N 2 i 1 1 + d(C i , C j ) d+1 4 ≤ C sup 1 n 2 n2 i=1 1 n 2 n2 j=1 ∆ 2 1 1 + d(C i , C j ) d+1 (Lemma A.27) ≤ C sup ∆ 2 n 2 max i=1,...,n2 n2 j=1 1 1 + d(C i , C j ) d+1 ≤ C sup ∆ 2 n 2 (Lemma A.
28, and because we will set ∆ → n→∞ ∞).

For T 4 in (11), using Cauchy-Schwarz and Lemma A.20,

T 4 = 1 n 2 n2 i=1 1 n 2 n2 j=1 1 ∆ 2 Ci dt i Cj dt j E r t θ (t i , X i ) R-1 i,θ R 0 (X i , X j ) R-1 j,θ r θ (t j , X j ) 2 ≤ C sup 1 n 2 n2 i=1 1 n 2 n2 j=1 1 ∆ 2 Ci dt i Cj dt j E [||r t θ (t i , X i )|| 4 ] E ||R 0 (X i , X j ) R-1 j,θ r θ (t j , X j )|| 4 . (14)
Using Condition 3.1, Lemma A.20 and Lemma A.29, we obtain

||R 0 (X i , X j ) R-1 j,θ r θ (t j , X j )|| 2 ≤ C sup N i N j 1 1 + d(C i , C j ) d+1 2 || R-1 j,θ r θ (t j , X j )|| 2 ≤ C sup N i N 2 j 1 1 + d(C i , C j ) d+1 2 .
Hence, going back to (14), .27 and A.28).

T 4 ≤ C sup 1 n 2 n2 i=1 1 n 2 n2 j=1 1 ∆ 2 Ci dt i Cj dt j E [N 2 i ] 1 1 + d(C i , C j ) d+1 2 E N 2 i N 4 j ≤ C sup 1 n 2 n2 i=1 1 n 2 n2 j=1 1 ∆ 2 Ci dt i Cj dt j 1 1 + d(C i , C j ) d+1 2 E [N 2 i ] E [N 4 i ] E N 8 j ≤ C sup ∆ 4 n 2 (Lemmas A
Hence, from (12), ( 13) and (15), we can set ∆ = n 1/6 to complete the proof.

Lemma A.12. For any fixed θ ∈ Θ,

E (|E (E n,θ ) -E [E n,θ |X]|)
goes to 0 as n → ∞.

Proof of Lemma A.12. Fix θ ∈ Θ. Because of Lemma A.10 and of |E(.)| ≤ E(|.|), it is sufficient to show that there exists a sequence ∆ → +∞ so that the lemma holds with E n,θ replaced by Ẽn,θ . Then, because of (E(.))

2 ≤ E (.) 2 , it is sufficient to show var E Ẽn,θ X → n→∞ 0.

Let C 1 , ..., C n2 be as in Definition A.7 and let ỹθ (t) be as in the proof of Lemma A.10. Define, for k = 1, ..., n 2 ,

g k (X) = 1 ∆ C k dt k E Y (t k ) -ỹθ (t k ) 2 X .
[Note that, following the discussion in Definition A.7, we have

g k (x) = K 0 (0) if N k = 0.] Then E Ẽn,θ |X = (1/n 2 ) n2 k=1 g k (X)
. Following the notation of Lemma A.11 we have, Let us also consider a Gaussian variable n+1 with mean 0 and variance δ 0 . Consider that X n+1 and n+1 are independent and independent of X, Y and . With the same argument as in the proof of Lemma A.10, we have

g k (X) = 1 ∆ C k dt k E Y (t k ) -r t θ (t k , X k ) R-1 k,θ y k 2 X ≤ 2 1 ∆ C k dt k K 0 (0) + r t θ (t k , X k ) R-1 k,θ Rk,0 R-1 k,θ r θ (t k , X k ) ≤ C sup + 2 1 ∆ C k dt k ||r t θ (t k , X k )|||| Rk,0 R-1 k,θ r θ (t k , X k )|| (Lemma A.20) ≤ C sup + C sup 1 ∆ C k dt k N k N k || R-1 k,θ
E (E n,θ ) = E [Y (X n+1 ) -ŷθ (X n+1 )] 2 ,
where we remind that ŷθ (X n+1 ) = r t θ (X n+1 )R -1 θ y, with r θ (X n+1 ) = (K(X 1 , X n+1 ), ..., K(X n , X n+1 )) t . Now, by symmetry of the roles of X 1 , ..., X n+1 and 1 , ... n+1 , we have

E (CV θ ) = E [Y (X n+1 ) -ŷn-1,θ (X n+1 )]
2 + δ 0 , with ŷn-1,θ (X n+1 ) = r t n-1,θ R-1 n-1,θ y, with r n-1,θ = (K(X 1 , X n+1 ), ..., K(X n-1 , X n+1 ), 0) t and Rn-1,θ = (K θ (X i , X j )) i,j=1,...,(n-1) + δ θ I n-1 0 0 1 .

Hence, using Cauchy-Schwarz

|E(CV θ ) -E(E n,θ ) -δ 0 | ≤ E r t θ (X n+1 )R -1 θ y -r t n-1,θ R-1 n-1,θ y 2 E r t θ (X n+1 )R -1 θ y + r t n-1,θ R-1 n-1,θ y -2Y (X n+1 ) 2 . ( 16
)
The second term in ( 16) is shown to be bounded with techniques similar to but simpler than in the proof of Lemma A.4. The first term in ( 16) is shown to go to zero with techniques similar to but simpler than in the proof of Lemma A.8. The five terms in the right-hand size of the above equation go to 0 in probability. Indeed, for fixed θ, the functions of θ go to 0 in probability because of Lemmas A.6, A.9, A.11, A.12 and A.13. The convergence of the supremums over θ to 0 is then a consequence of the fact that Θ is compact and of Lemma A.4 and corollaries A.5 and A.14. Finally, since θCV minimizes CV θ + δ 0 , we conclude with, for any θ ∈ Θ Lemma A.25. Use the notation n 2 , ∆ and C 1 , ..., C n2 of Definition A.7. Let, for i = 1, ..., n 2 , X i 1 , ..., X i Ni be the N i components of X that are in C i (so that the order of their indexes in X is preserved). Then i) For i = 1, ..., n 2 , N i follows a binomial B(n, 1/n 2 ) distribution. For any i, j = 1, ..., n 2 ; i = j, conditionally to N i = k i , N j follows a binomial B(n -k i , 1/(n 2 -1)) distribution.

Corollary

ii) Conditionally to N i = k i , X i 1 , ..., X i ki are independent and uniformly distributed on C i .

iii) For 1 ≤ i = j ≤ n 2 , conditionally to N i = k i , N j = k j , the sets of random variables (X i 1 , ..., X i ki ) and (X j 1 , ..., X j kj ) are independent, and their components are independent and uniformly distributed on C i and C j respectively.

Consider n 2 real-valued functions f 1 , ..., f n2 of X that can be written f i (X) = f (N i , X i 1 , ..., X i Ni ), and so that, for any t ∈ R d , x 1 , ..., x N ∈ R d , f (N, x 1 + t, ..., x N + t) = f (N, x 1 , ..., x N ). Then iv) The variables f 1 (X), ..., f n2 (X) have the same distribution. The couples (f i (X), f j (X)), for 1 ≤ i = j ≤ n 2 , have the same distribution.

Lemma A.26. Use the notation of Lemma A.25, and consider n 2 functions f 1 , ..., f n2 that satisfy the conditions of Lemma A.25. Assume that there exist fixed even natural numbers q, l and a finite constant C sup (independent of n and X) so that E f 2 i (X)|N i = k ≤ C sup (1 + k q + k q+l /∆ l ). Then, if ∆ → n→∞ +∞ and ∆ = O(n 1/(2q+5) ),

var 1 n 2 n2 i=1 f i (X) → n→∞ 0.
Lemma A.27. Let N follow the binomial distribution B(n, 1/n 2 ), with n/n 2 = ∆ → n→∞ +∞. Then, for any k ∈ N, there exists a finite constant C sup , independent of n, so that

E N k ≤ C sup ∆ k .
Lemma A.28. Let n 2 , ∆ and C 1 , ..., C n2 be as in Definition A.7. Assume that ∆ is lower bounded, as a function of n. Then, there exists a finite constant C sup so that for any n, i ∈ {1, ..., n 2 },

n2 j=1 1 1 + d(C i , C j ) d+1 ≤ C sup .
Lemma A.29. Let A be a real m 1 × m 2 matrix and b be a m 2 -dimensional real column vector. Then

||Ab|| 2 ≤ m 1 m 2 max i,j A 2 i,j ||b|| 2 .

Figure 1 :

 1 Figure1: Simulation of N = 1000 independent realizations of n = 100 iid observation points with uniform distribution on [0, n], of the Gaussian process Y on [0, n] with Matérn (σ 2 0 = 1, 0 = 3, ν 0 = 10) covariance function, and of the corresponding observation vector with noise variance δ 0 = 0.25 2 . For each simulation, ν 0 is known, the noise variance is fixed to δ 1 = δ 0 (well-specified case) or δ 1 = 0.1 2 = δ 0 (misspecified case), σ 2 and are estimated by ML and CV and the corresponding error criteria D n,σ 2 , ˆ (normalized Kullback-Leibler divergence) and E n,σ 2 , ˆ (integrated square prediction error) are computed. The histograms of the N estimates of and of the N corresponding values of the error criteria are reported for ML and CV and in the well-specified and misspecified cases. In the well-specified case, the estimates are on average reasonably close to the true values, the error criteria are reasonably small and ML performs better than CV in all aspects. In the misspecified case, the ML and CV estimates of the correlation lengths are significantly different, ML performs better than CV for D n,σ 2 , ˆ and CV performs better than ML for E n,σ 2 , ˆ .

Figure 2 :

 2 Figure2: Same settings as for Figure1but with N = 200 and n = 500. The conclusions for the comparison of ML and CV are the same as for Figure1. The variances of the estimates ˆ and of the error criteria are smaller than in Figure1and the histograms of the estimates are approximately symmetric and unimodal. Between Figure1and Figure2, all error criteria keep the same order of magnitude except D n,σ 2 , ˆ in the well-specified case which decreases and becomes very small for ML and CV (see also Table1). This is because in the well-specified case D n,σ 2 0 , 0 = 0.

  r θ (t k , X k )|| (Conditions 3.1 and 3.2 and Lemma A.29)≤ C sup 1 + N 2 k (Lemma A.20). Hence E g 2 k (X) N k = N ≤ C sup (1 + N 4), so that we can complete the proof with Lemma A.26.Lemma A.13. Consider a fixed θ ∈ Θ. ThenE(CV θ ) -E(E n,θ ) -δ 0 goes to 0 as n → ∞.Proof of Lemma A.13. Let us consider a random observation point X n+1 with uniform distribution on [0, n 1/d ] d .

  considering t as a new random observation point X n+1 as in the proof of Lemma A.10, we show the first bound of the lemma as in the proof of Lemma A.4, the only difference being that there are n + 1 observation points instead of n. The second and third bounds are proved as in the proof of corollary A.5.Proof of Theorem 3.4. We havesup θ∈Θ |CV θ -δ 0 -E n,θ | ≤ sup θ∈Θ |CV θ -E(CV θ |X)| + sup θ∈Θ |E(CV θ |X) -E(CV θ )| + sup θ∈Θ |E(CV θ ) -δ 0 -E(E n,θ )| + sup θ∈Θ |E(E n,θ ) -E(E n,θ |X)| + sup θ∈Θ |E(E n,θ |X) -E n,θ | .

E

  n, θCV -E n,θ ≤ CV θCV -CV θ + 2 sup θ∈Θ |CV θ -δ 0 -E n,θ | ≤ 2 sup θ∈Θ |CV θ -δ 0 -E n,θ | .Hence sup θ∈Θ E n, θCV -E n,θ = o p (1).

Table 1 :

 1 For the settings of Figures1 and 2, the averages and standard-deviations of ˆ , D n,σ 2 , ˆ and E n,σ 2 , ˆ are reported.

	n	Specification Estimation Average Standard deviation Average	Average
				of ˆ	of ˆ	of E n,σ 2 , ˆ of D n,σ 2 , ˆ
		Well-specified	ML	3.031	0.370	0.073	0.023
	100	Well-specified Misspecified	CV ML	3.447 1.090	1.159 0.553	0.085 0.255	0.222 1.025
		Misspecified	CV	3.525	1.324	0.087	3.563
		Well-specified	ML	3.003	0.177	0.070	0.005
	500	Well-specified Misspecified	CV ML	3.084 0.975	0.399 0.269	0.071 0.247	0.037 0.972
		Misspecified	CV	3.089	0.437	0.071	3.442

  j and ||A|| for the largest singular value of A. Both |.| and ||.|| are norms and ||.|| is also a matrix norm.

  A.14. For any i = 1, ..., p, For fixed n we can exchange derivative and integration, so we obtain

	E sup θ∈Θ	∂ ∂θ i	E n,θ , E sup θ∈Θ	∂ ∂θ i	E ( E n,θ | X)	and sup θ∈Θ	∂ ∂θ i	E (E n,θ )
	are bounded w.r.t n.
	Proof of corollary A.14. We have
	E sup θ∈Θ		∂ ∂θ i	E n,θ	= E sup θ∈Θ	∂ ∂θ i	1 n [0,n 1/d ] d	[Y (t) -ŷθ (t)] 2 dt .
	E sup θ∈Θ	∂ ∂θ i	E n,θ	= E sup θ∈Θ	1 n [0,n 1/d ] d	∂ ∂θ i	[Y (t) -ŷθ (t)] 2 dt
						≤ E	1 n [0,n 1/d ] d	sup

θ∈Θ ∂ ∂θ i [Y (t) -ŷθ (t)]

Acknowledgements

The author thanks Josselin Garnier and Benedikt Pötscher for constructive discussions.

Let us use the notation C 1 , ..., C n2 of Definition A.7. Let, for t ∈ R d and v = (v 1 , ..., v m ) ∈ (R d ) m , r θ (t, v) = (K θ (t, v 1 ), ..., K θ (t, v m )) t . We define r 0 (t, v) similarly. Let y i , Ri,θ and Ri,0 be as in Definition A.7. Let for i = j, R 0 (X i , X j ) = K 0 ((X i ) k , (X j ) l ) k=1,...,Ni;l=1,...,Nj . Let R 0 (X i , X i ) = K 0 ((X i ) k , (X i ) l ) k,l=1,...,Ni + δ 0 I Ni . Then,

Hence, using the relation cov A 2 , B 2 = 2 (cov(A, B)) 2 , for two centered Gaussian variables A and B, we obtain

where T 1 , T 2 , T 3 , T 4 are defined and treated below, and with T 2 = T 3 by symmetry.

For T 1 ,

For T 2 , using Cauchy-Schwarz and Lemma A.20,

A.4 Technical results

The following technical results are proved in the supplementary material.

Lemma A.15. Consider a fixed number n of observation points. Consider a function f θ (X, y) that is p times continuously differentiable w.r.t θ for any X, y and so that, for i 1 + ...

has finite mean value w.r.t X and y. Then, there exists a constant C sup (depending only of Θ) so that

Lemma A.16. There exists a finite constant C sup so that, for any a, b ∈ R d ,

Lemma A.17. Let 0 < C inf ≤ C sup < ∞ be fixed independently of n. Let s n be a function of n so that s n ∈ N * and C inf n ≤ s n ≤ C sup n. Consider s n observation points X1 , ..., Xsn , independent and uniformly distributed on [0, n 1/d ] d . Let A 1 , ..., A k be k sequences of s n × s n random matrices so that, for l = 1, ..., k, (A l ) i,j depends only on Xi and Xj and satisfies

Lemma A.18. The supremum over n, θ and X of the eigenvalues of

Lemma A.19. Lemma A.18 also holds when K θ is replaced by Kθ of Definition A.7.

Lemma A.20. Lemma A.18 also holds when R θ is replaced by Rk,θ of Definition A.7.

Lemma A.21. Let k ∈ N. Let A 1,θ , ..., A k,θ be k sequences of symmetric random matrices (functions of X and θ) so that, for any m ∈ N, a 1 , ..., a m ∈ {1, ..., k}, sup θ∈Θ E X |A a1,θ ...A am,θ | 2 is bounded (w.r.t n). Let B 1,θ , ..., B k+1,θ be k + 1 sequences of random symmetric non-negative matrices (functions of X and θ) so that sup θ ||B 1,θ ||, ..., sup θ ||B k+1,θ || are bounded (w.r.t n and X). Then

Lemma A.22. Consider a fixed θ ∈ Θ. With the notation of Definition A.7, we have, when

Then, there exists a finite constant C sup so that, for any n, for any

Lemma A.24. Use the notation n 2 , ∆, C(t), f (T ) and D ∆ (x 1 , x 2 ) of Definition A.7 and Lemma A.23. Then, when n 2 = o(n),

Supplementary material

In the supplementary material, we give the proof of the lemmas stated in Section A.4.