Confidence bands for HorvitzThompson estimators
using sampled noisy functional data
Abstract
When collections of functional data are too large to be exhaustively observed, survey sampling techniques provide an effective way to estimate global quantities such as the population mean function. Assuming functional data are collected from a finite population according to a probabilistic sampling scheme, with the measurements being discrete in time and noisy, we propose to first smooth the sampled trajectories with local polynomials and then estimate the mean function with a HorvitzThompson estimator. Under mild conditions on the population size, observation times, regularity of the trajectories, sampling scheme, and smoothing bandwidth, we prove a Central Limit Theorem in the space of continuous functions. We also establish the uniform consistency of a covariance function estimator and apply the former results to build confidence bands for the mean function. The bands attain nominal coverage and are obtained through Gaussian process simulations conditional on the estimated covariance function. To select the bandwidth, we propose a crossvalidation method that accounts for the sampling weights. A simulation study assesses the performance of our approach and highlights the influence of the sampling scheme and bandwidth choice.
Keywords : CLT, functional data, local polynomial smoothing, maximal inequalities, space of continuous functions, suprema of Gaussian processes, survey sampling, weighted crossvalidation.
1 Introduction
The recent development of automated sensors has given access to very large collections of signals sampled at fine time scales. However, exhaustive transmission, storage, and analysis of such massive functional data may incur very large investments. In this context, when the goal is to assess a global indicator like the mean temporal signal, survey sampling techniques are appealing solutions as they offer a good tradeoff between statistical accuracy and global cost of the analysis. In particular they are competitive with signal compression techniques (Chiky and Hébrail, 2008). The previous facts provide some explanation why, although survey sampling and functional data analysis have been longestablished statistical fields, motivation for studying them jointly only recently emerged in the literature. In this regard Cardot et al. (2010a) examine the theoretical properties of functional principal components analysis (FPCA) in the survey sampling framework. Cardot et al. (2010b) harness FPCA for modelassisted estimation by relating the unobserved principal component scores to available auxiliary information. Focusing on sampling schemes, Cardot and Josserand (2011) estimate the mean electricity consumption curve in a population of about 19,000 customers whose electricity meters were read every 30 minutes during one week. Assuming exact measurements, they first perform a linear interpolation of the discretized signals and then consider a functional version of the HorvitzThompson estimator. For a fixed sample size, they show that estimation can be greatly improved by utilizing stratified sampling over simple random sampling and they extend the Neyman optimal allocation rule (see e.g. Särndal et al. 1992) to the functional setup. Note however that the finitesample and asymptotic properties of their estimator rely heavily on the assumption of errorfree measurements, which is not always realistic in practice. The first contribution of the present work is to generalize the framework of Cardot and Josserand (2011) to noisy functional data. Assuming curve data are observed with errors that may be correlated over time, we replace the interpolation step in their procedure by a smoothing step based on local polynomials. As opposed to interpolation, smoothing can effectively reduce the noise level in the data, which improves estimation accuracy. We establish a functional CLT for the mean function estimator based on the smoothed data and prove the uniform consistency of a related covariance estimator. These results have important applications to the simultaneous inference of the mean function.
In relation to mean function estimation, a key statistical task is to build confidence regions. There exists a vast and still active literature on confidence bands in nonparametric regression. See e.g. Sun and Loader (1994), Eubank and Speckman (1993), Claeskens and van Keilegom (2003), Krivobokova et al. (2010), and the references therein. When data are functional the literature is much less abundant. One possible approach is to obtain confidence balls for the mean function in a space. Mas (2007) exploits this idea in a goodnessoffit test based on the functional sample mean and regularized inverse covariance operator. Using adaptive projection estimators, Bunea et al. (2011) build conservative confidence regions for the mean of a Gaussian process. Another approach consists in deriving results in a space of continuous functions equipped with the supremum norm. This allows for the construction of confidence bands that can easily be visualized and interpreted, as opposed to confidence balls. This approach is adopted for example by Faraway (1997) to build bootstrap bands in a varyingcoefficients model, by Cuevas et al. (2006) to derive bootstrap bands for functional location parameters, by Degras (2009, 2011) to obtain normal and bootstrap bands using noisy functional data, and by Cardot and Josserand (2011) in the context of a finite population. In the latter work, the strategy was to first establish a CLT in the space and then derive confidence bands based on a simple but rough approximation to the supremum of a Gaussian process (Landau and Shepp, 1970). Unfortunately, the associated bands depend on the datagenerating process only through its variance structure and not its correlation structure, which may cause the empirical coverage to differ from the nominal level. The second innovation of our paper is to propose confidence bands that are easy to implement and attain nominal coverage in the survey sampling/finite population setting. To do so we use Gaussian process simulations as in Cuevas et al. (2006) or Degras (2011). This procedure can be thought as a parametric bootstrap, where the parameter to be estimated, the covariance function, is lying in an infinite dimensional functional space. Our contribution is to provide the theoretical underpinning of the construction method, thereby guaranteeing that nominal coverage is attained asymptotically. The theory we derive involves maximal inequalities, random entropy numbers, and large covariance matrix theory.
Finally, the implementation of the mean function estimator developed in this paper requires the selection of a bandwidth in the data smoothing step. Objective, datadriven bandwidth selection methods are desirable for this purpose. As explained by Opsomer and Miller (2005), bandwidth selection in the survey estimation context poses specific problems (in particular, the necessity to take the sampling design into account) that make usual crossvalidation or mean square error optimization methods inadequate. In view of the modelassisted survey estimation of a population total, these authors propose a crossvalidation method that aims at minimizing the variance of the estimator, the bias component being negligible in their setting. In our functional and designbased framework, the bias is however no longer negligible. We therefore devise a novel crossvalidation criterion based on weighted least squares, with weights proportional to the sampling weights. For the particular case of simple random sampling without replacement, this criterion reduces to the cross validation technique of Rice and Silverman (1991), whose asymptotic properties has been studied by Hart and Wehrly (1993).
The paper is organized as follows. We fix notations and define our estimators in section 2. In section 3, we introduce our asymptotic framework based on superpopulation models (see Isaki and Fuller, 1982), establish a CLT for the mean function estimator in the space of continuous functions, and show the uniform consistency of a covariance estimator. Based on these results, we propose a simple and effective method for building simultaneous confidence bands. In section 4, a weighted crossvalidation procedure is proposed for selecting the bandwidth and simulations are performed to compare different sampling schemes and bandwidth choices. Our estimation methodology is seen to compare favorably with other methods and to achieve nearly optimal performances. The paper ends with a short discussion on topics for future research. Proofs are gathered in an Appendix.
2 Notations and estimators
Consider a finite population of size and suppose that to each unit corresponds a real function on with We assume that each trajectory belongs to the space of continuous functions Our target is the mean trajectory defined as follows:
(1) 
We consider a random sample drawn from without replacement according to a fixedsize sampling design where is the probability of drawing the sample The size of is nonrandom and we suppose that the first and second order inclusion probabilities satisfy

for all

for all
so that each unit and each pair of units can be drawn with a non null probability from the population. Note that for simplicity of notation the subscript has been omitted. Also, by convention, we write for all .
Assume that noisy measurements of the sampled curves are available at fixed discretization points For all units , we observe
(2) 
where the measurement errors are centered random variables that are independent across the index (units) but not necessarily across (possible temporal dependence). It is also assumed that the random sample is independent of the noise and the trajectories are deterministic.
Our goal is to estimate as accurately as possible and to build asymptotic confidence bands, as in Degras (2011) and Cardot and Josserand (2011). For this, we must have a uniformly consistent estimator of its covariance function.
2.1 Linear smoothers and the HorvitzThompson estimator
For each (potentially observed) unit , we aim at recovering the curve by smoothing the corresponding discretized trajectory with a linear smoother (e.g. spline, kernel, or local polynomial):
(3) 
Note that the reconstruction can only be performed for the observed units .
Here we use local linear smoothers (see e.g. Fan and Gijbels (1997)) because of their wide popularity, good statistical properties, and mathematical convenience. The weight functions can be expressed as
(4) 
where is a kernel function, is a bandwidth, and
(5) 
We suppose that the kernel is nonnegative, has compact support, satisfies and for some finite constant and for all .
The classical HorvitzThompson estimator of the mean curve is
(6)  
where is the sample membership indicator ( if and otherwise). It holds that and
2.2 Covariance estimation
The covariance function of can be written as
(7) 
for all , where
(8) 
with
(9) 
A natural estimator of is given by
(10) 
It is unbiased and its uniform mean square consistency is established in Section 3.2.
3 Asymptotic theory
We consider the superpopulation framework introduced by Isaki and Fuller (1982) and discussed in detail by Fuller (2009). Specifically, we study the behaviour of the estimators and as population increases to infinity with . Recall that the sample size , inclusion probabilities and , and grid size all depend on . In what follows we use the notations and for finite, positive constants whose value may vary from place to place. The following assumptions are needed for our asymptotic study.

(Sampling design) for all () and . and

(Trajectories) and for all and , where is a finite constant.

(Growth rates) for all and . as

(Measurement errors) The random vectors are i.i.d. and follow the multivariate normal distribution with mean zero and covariance matrix . The largest eigenvalue of the covariance matrix satisfies for all .
Assumption (A1) deals with the properties of the sampling design. It states that the sample size must be at least a positive fraction of the population size, that the one and twofold inclusion probabilities must be larger than a positive number, and that the twofold inclusion probabilities should not be too far from independence. The latter is fulfilled for example for stratified sampling with sampling without replacement within each stratum (Robinson and Särndal (1983)) and is discussed in details in Hàjek (1981) for rejective sampling and other unequal probability sampling designs. Assumption (A2) imposes Hölder continuity on the trajectories, a mild regularity condition. Assumption (A3) states that the design points have a quasiuniform repartition (this holds in particular for equidistant designs and designs generated by a regular density function) and that the grid size is essentially negligible compared to the population size (for example if for some ). In fact the results of this paper also hold if stays bounded away from zero and infinity as (see Section 5). Finally (A4) imposes joint normality, short range temporal dependence, and bounded variance for the measurement errors . It is trivially satisfied if the are independent with variances . It is also verified if the arise from a discrete time Gaussian process with short term temporal correlation such as ARMA or stationary mixing processes. Note that the Gaussian assumption is not central to our derivations: it can be weakened and replaced by moment conditions on the error distributions at the expense of much more complicated proofs.
3.1 Limit distribution of the HorvitzThompson estimator
We now derive the asymptotic distribution of our estimator in order to build asymptotic confidence bands. Obtaining the asymptotic normality of estimators in survey sampling is a technical and difficult issue even for simple quantities such as means or totals of real numbers. Although confidence intervals are commonly used in the survey sampling community, the Central Limit Theorem (CLT) has only been checked rigorously, as far as we know, for a few sampling designs. Erdös and Rényi (1959) and Hàjek (1960) proved that the HorvitzThompson estimator is asymptotically Gaussian for simple random sampling without replacement. The CLT for rejective sampling is shown by Hàjek (1964) whereas the CLT for other proportional to size sampling designs is studied by Berger (1998). Recently these results were extended for some particular cases of twophase sampling designs (Chen and Rao, 2007). Let us assume that the HorvitzThompson estimator satisfies a CLT for realvalued quantities.

(Univariate CLT) For any fixed , it holds that
as , where stands for convergence in distribution.
We recall here the definition of the weak convergence in equipped with the supremum norm (e.g. van der Vaart and Wellner, 2000). A sequence of random elements of is said to converge weakly to a limit in if as for all bounded, uniformly continuous functionals on .
To establish the limit distribution of in , we need to assume the existence of a limit covariance function
In the following theorem we state the asymptotic normality of the estimator in the space equipped with the sup norm.
Theorem 1.
Assume (A1)–(A5) and that and as . Then
in , where is a Gaussian process with mean zero and covariance function .
Theorem 1 provides a convenient way to infer the local features of . It is applied in Section 3.3 to the construction of simultaneous confidence bands, but it can also be used for a variety of statistical tests based on supremum norms (see Degras, 2011).
Observe that the conditions on the bandwidth and design size are not very constraining. Suppose for example that and for some . Then and satisfy the conditions of Theorem 1 as soon as . Thus, for more regular trajectories, i.e. larger the bandwidth can be chosen with more flexibility.
The proof of Theorem 1 is similar in spirit to that of Theorem 1 in Degras (2011) and Proposition 3 in Cardot and Josserand (2010). Essentially, it breaks down into: (i) controlling uniformly on the bias of , (ii) establishing the functional asymptotic normality of the local linear smoother applied to the sampled curves , and (iii) controlling uniformly on (in probability) the local linear smoother applied to the errors . Part (i) is easily handled with standard results on approximation properties of local polynomial estimators (see e.g. Tsybakov (2009)). Part (ii) mainly consists in proving an asymptotic tightness property, which entails the computation of entropy numbers and the use of maximal inequalities (van der Vaart and Wellner, 2000). Part (iii) requires first to show the finitedimensional convergence of the smoothed error process to zero and then to establish its tightness with similar arguments as in part (ii).
3.2 Uniform consistency of the covariance estimator
We first note that under (A1)–(A4), by the approximation properties of local linear smoothers, converges uniformly to on as and . Hence the consistency of can be stated with respect to instead of . In alignment with the related Proposition 2 in Cardot and Josserand (2011) and Theorem 3 in Breidt and Opsomer (2000), we need to make some assumption on the twofold inclusion probabilities of the sampling design :

where is the set of all quadruples in with distinct elements.
This assumption is discussed in detail in Breidt and Opsomer (2000) and is fulfilled for example for stratified sampling.
Theorem 2.
Assume (A1)–(A4), (A6), and that and for some as . Then
where the expectation is jointly with respect to the design and the multivariate normal model.
Note the additional condition on the bandwidth in Theorem 2. If we suppose, as in the remark in Section 3.1, that and for some , then condition as is fulfilled with e.g. .
3.3 Confidence bands
In this section we build confidence bands for of the form
(11) 
where is a suitable number and . More precisely, given a confidence level we seek that approximately satisfies
(12) 
where is a Gaussian process with mean zero and covariance function , and where . Exact bounds for the supremum of Gaussian processes have been derived for only a few particular cases (Adler and Taylor, 2007, Chapter 4). Computing accurate and as explicit as possible bounds in a general setting is a difficult issue and would require additional strong conditions such as stationarity which have no reason to be fulfilled in our setting.
In view of Theorems 12 and Slutski’s Theorem, the bands defined in (11) with chosen as in (12) will have approximate coverage level . The following result provides a simulationbased method to compute .
Theorem 3.
Assume (A1)–(A6) and for some as . Let be a Gaussian process with mean zero and covariance function . Let be a sequence of processes such that for each , conditionally on , is Gaussian with mean zero and covariance defined in (10). Then for all , as , the following convergence holds in probability:
Theorem 3 is derived by showing the weak convergence of to in , which stems from Theorem 2 and the Gaussian nature of the processes . As in the first two theorems, maximal inequalities are used to obtain the above weak convergence. The practical importance of Theorem 3 is that it allows to estimate the number in (12) via simulation: (with the previous notations), conditionally on , one can simulate a large number of sample paths of the Gaussian process and compute their supremum norms. One then obtains a precise approximation to the distribution of , and it suffices to set as the quantile of order of this distribution:
(13) 
4 A simulation study
In this section, we evaluate the performances of the mean curve estimator as well as the coverage and the width of the confidence bands for different bandwidth selection criteria and different levels of noise. The simulations are conducted in the R environment.
4.1 Simulated data and sampling designs
We have generated a population of curves discretized at and equidistant instants of time in The curves of the population are generated so that they have approximately the same distribution as the electricity consumption curves analyzed in Cardot and Josserand (2011) and each individual curve for is simulated as follows
(14) 
where the mean function is drawn in Figure 2 and the random variables are independent realizations of a centered Gaussian random variable with variance . The three basis function and are orthonormal functions which represent the main mode of variation of the signals, they are represented in Figure 1. Thus, the covariance function of the population is simply
(15) 
To select the samples, we have considered two probabilistic selection procedures, with fixed sample size,

Simple random sampling without replacement (SRSWOR).

Stratified sampling with SRSWOR in all strata. The population is divided into a fixed number of strata built by considering the quantiles and of the total consumption for all units . For example, the first strata contains all the units such that and thus its size is half of the population size The sample size in stratum is determined by a Neymanlike allocation, as suggested in Cardot and Josserand (2011), in order to get a HorvitzThompson estimator of the mean trajectory whose variance is as small as possible. The sizes of the different strata, which are optimal according to this mean variance criterion, are reported in Table 1.
Stratum number  1  2  3  4  5 

Stratum size  10000  4000  3000  2000  1000 
Allocation  655  132  98  68  47 
We suppose we observe, for each unit in the sample the discretized trajectories, at equispaced points,
(16) 
The parameter controls the noise level compared to the true signal. We consider two different situations for the noise components

Heteroscedasticity. The are independent random variables whose variances are proportional to the population variances at time .

Temporal dependence. The are stationary AR(3) processes with Gaussian innovations generated as follows
The are i.i.d. and is such that
As an illustrative example, a sample of noisy discretized curves are plotted in Figure 2 with heteroscedastic noise components and in Figure 3 for correlated noise. It should be noted that the observed trajectories corrupted by the correlated noise are much smoother than the trajectories corrupted by the heteroscedastic noise. The empirical standard deviation in the population, for these two different type of noise are drawn in Figure 4.
4.2 Weighted crossvalidation for bandwidth selection
Assuming we can access the exact trajectories (which is the case in simulations) we consider the oracletype estimator
(17) 
which will be a benchmark in our numerical study. We compare different interpolation and smoothing strategies for estimating the :

Linear interpolation of the as in Cardot and Josserand (2011).

Local linear smoothing of the with bandwidth as in (3).
The crucial parameter here is . To evaluate the interest of smoothing and the performances of datadriven bandwidth selection criteria, we consider an error measure that compares the oracle to any estimator based on the noisy data :
(18) 
Considering the estimator defined in (6), we denote by the bandwidth that minimizes (18). The mean estimator built with bandwidth is called smooth oracle estimator.
When , as in SRSWOR and stratified sampling, it can be easily checked that is the minimum argument of the weighted least squares functional
(19) 
with respect to where the weights are Then, a simple and natural way to select bandwidth is to consider the following designbased cross validation
(20) 
where
with new weights A heuristic justification for this approach is that, given we have for and Thus,
and, up to which does not depend on , the minimum value of the expected cross validation criterion should be attained for estimators which are not too far from
This weighted cross validation criterion is simpler than the cross validation criteria based on the estimated variance proposed in Opsomer and Miller (2005). Indeed, in our case, the bias may be non negligible and focusing only on the variance part of the error leads to too large selected values for the bandwidth. Furthermore, Opsomer and Miller (2005) suggested to consider weights defined as follows For SRSWOR, since one has so that the weighted cross validation criterion defined in (20) is exactly the cross validation criterion introduced by Rice and Silverman (1991) in the independent case. We denote in the following by the bandwidth value minimizing this criterion.
For stratified sampling, a better approximation which keeps the designbased properties of the estimator can be obtained by taking into account the sampling rates in the different strata. Assume the population is partitioned in strata of respective sizes and we sample observations in each by SRSWOR. If we have Thus, we take for all the units in stratum and scale the weights for all the units of the sample that do not belong to stratum We denote by the bandwidth value minimizing (20).
4.3 Estimation errors and confidence bands
We draw samples in the population of curves and compare the different estimators of Section 4.2 with the loss criterion
(21) 
for different values of and in (16). For comparison, the quadratic approximation error for function by its average value, is
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  17.65  3.08  8.73  23.50  4.22  1.44  2.79  5.59 
17.65  3.07  8.71  23.51  6.49  3.61  5.36  8.03  
17.65  3.07  8.71  23.51  4.22  1.45  2.78  5.56  
17.65  3.07  8.72  23.50  4.22  1.45  2.78  5.57  
17.60  3.01  8.70  23.36  4.17  1.38  2.76  5.55  
25%  lin  17.69  3.94  8.99  21.52  5.26  2.63  4.15  6.54 
17.53  3.83  8.76  21.53  6.98  4.29  5.83  8.47  
17.53  3.83  8.76  21.53  5.02  2.39  3.89  6.33  
17.52  3.81  8.78  21.52  5.01  2.37  3.88  6.27  
16.58  2.85  7.87  20.01  4.07  1.46  2.94  5.28 
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  18.03  3.39  9.24  23.27  4.05  1.45  2.86  5.35 
18.02  3.38  9.26  23.34  6.09  3.24  4.87  7.56  
18.02  3.38  9.26  23.34  4.05  1.45  2.82  5.40  
18.02  3.38  9.27  23.32  4.04  1.43  2.83  5.39  
17.98  3.35  9.20  23.17  4.00  1.39  2.81  5.29  
25%  lin  18.16  3.89  9.43  22.86  5.25  2.85  4.24  6.57 
17.55  3.30  8.89  22.09  6.45  3.77  5.37  8.11  
17.55  3.30  8.89  22.09  4.57  2.12  3.49  5.81  
17.55  3.28  8.89  22.09  4.56  2.11  3.48  5.81  
17.04  2.75  8.38  21.87  4.04  1.60  3.02  5.31 
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  16.23  3.05  8.67  20.86  4.08  1.40  2.88  5.44 
16.24  3.07  8.66  20.88  5.90  2.99  4.70  7.33  
16.24  3.07  8.66  20.88  4.10  1.38  2.90  5.47  
16.24  3.06  8.65  20.88  4.10  1.38  2.90  5.46  
16.19  3.01  8.69  20.86  4.04  1.34  2.82  5.36  
25%  lin  17.18  3.88  9.38  22.04  5.22  2.65  4.07  6.47 
17.13  3.84  9.28  22.02  6.76  3.98  5.76  8.32  
17.13  3.84  9.28  22.02  5.16  2.59  4.02  6.37  
17.12  3.81  9.25  22.02  5.15  2.59  4.01  6.37  
16.12  2.87  8.17  21.00  4.04  1.49  2.94  5.27 
The empirical mean as well as the first, second and third quartiles of the estimation error are given, when d=200, in Table 2 for the heteroscedastic noise case. Results for are presented for the heteroscedastic case in Table 3 and in Table 4 for the correlated case.
We first note that in all simulations, stratified sampling largely improves the estimation of the mean curve in comparison to SRSWOR. Also, linear interpolation performs nearly as well as the smooth oracle estimator for large samples, especially when the noise level is low (). As far as bandwidth selection is concerned, the usual cross validation criterion is not adapted to unequal probability sampling and tends to select too large bandwidth values. In particular it does not perform as well as linear interpolation for stratified sampling. On the other hand, our weighted crossvalidation method seems effective for selecting the bandwidth. It produces estimators that are very close to the oracle and that dominate the other estimators when the noise level is moderate or high ().
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  0.044  0.041  0.044  0.047  0.049  0.046  0.049  0.053 
0.044  0.041  0.044  0.048  2.520  2.083  2.852  3.032  
0.044  0.041  0.044  0.048  0.058  0.054  0.058  0.062  
0.044  0.041  0.044  0.047  0.049  0.045  0.049  0.052  
25%  lin  1.087  1.011  1.080  1.156  1.214  1.134  1.210  1.287 
0.905  0.837  0.901  0.970  3.155  2.638  3.260  3.602  
0.905  0.837  0.901  0.970  1.009  0.936  1.004  1.076  
0.898  0.830  0.894  0.962  0.990  0.919  0.988  1.055 
This is clearer when we look at criterion defined in (18), which only focuses on the part of the estimation error which is due to the noise. Results are presented in Table 5 for in the heteroscedastic case. For errors are given in Table 6 in the heteroscedastic case and in Table 7 for correlated noise. When the noise level is high, we observe a significant impact of the number of discretization points on the accuracy of the smoothed estimators. Our individual trajectories, which have roughly the same shape as load curves, are actually not very smooth so that smoothing approaches are only really interesting, compared to linear interpolation, when the number of discretization points is large enough. Finally, it also becomes clearer that a key parameter is the bandwidth value which has to be chosen with appropriate criteria that must take the sampling weights into account. When the noise level is low (), the error according to criterion is multiplied by at least 15 in stratified sampling.
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  0.044  0.042  0.044  0.047  0.049  0.047  0.049  0.051 
0.040  0.038  0.040  0.042  2.231  1.612  1.917  2.806  
0.040  0.038  0.040  0.042  0.052  0.049  0.052  0.055  
0.040  0.038  0.040  0.042  0.044  0.041  0.044  0.046  
25%  lin  1.089  1.030  1.087  1.142  1.219  1.155  1.212  1.280 
0.498  0.462  0.495  0.535  2.591  1.932  2.344  3.254  
0.498  0.462  0.495  0.535  0.552  0.509  0.549  0.594  
0.497  0.460  0.494  0.533  0.547  0.505  0.545  0.586 
SRSWOR  Stratified sampling  
Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  0.17  0.15  0.16  0.18  0.05  0.04  0.04  0.05 
0.17  0.15  0.16  0.18  1.94  1.53  1.59  2.90  
0.17  0.15  0.16  0.18  0.07  0.07  0.07  0.08  
0.17  0.15  0.16  0.18  0.07  0.07  0.07  0.08  
25%  lin  1.09  1.03  1.09  1.14  1.20  1.08  1.19  1.32 
0.50  0.46  0.50  0.53  2.83  2.19  2.57  3.67  
0.50  0.46  0.50  0.53  1.15  1.02  1.13  1.26  
0.49  0.46  0.49  0.53  1.13  1.01  1.12  1.25 
SRSWOR  Stratified sampling  

Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  97.2  10.91  10.74  10.90  11.07  98.1  5.95  5.87  5.95  6.02 
97.3  10.89  10.73  10.89  11.06  47.5  5.68  5.60  5.68  5.76  
97.3  10.89  10.73  10.89  11.06  97.5  5.92  5.84  5.91  6.00  
97.2  10.90  10.72  10.90  11.07  98.0  5.94  5.86  5.94  6.02  
97.3  10.54  10.36  10.54  10.70  98.2  5.59  5.51  5.60  5.67  
25%  lin  97.7  13.23  13.06  13.22  13.41  98.3  8.27  8.19  8.27  8.36 
97.2  12.66  12.49  12.65  12.83  64.7  6.70  6.60  6.69  6.79  
97.2  12.66  12.49  12.65  12.83  97.3  7.56  7.48  7.56  7.65  
97.3  12.70  12.50  12.70  12.87  97.5  7.68  7.58  7.68  7.79  
97.0  10.53  10.37  10.52  10.70  97.7  5.59  5.51  5.59  5.66 
We now examine in Table 8, Table 9 and Table 10 the empirical coverage and the width of the confidence bands, which are built as described in Section 3.3. For each sample, we estimate the covariance function and draw realizations of a centered Gaussian process with variance function in order to obtain a suitable coefficient with a confidence level of as explained in equation (13). The area of the confidence band is then The results highlight now the interest of considering smoothing strategies combined with the weighted cross validation bandwidth selection criterion (20). For stratified sampling, the use of the unweighted cross validation criterion leads to empirical coverage levels that are significantly below the nominal one. It also appears that linear interpolation, which does not intend to get rid of the noise, always gives larger confidence bands than the smoothed estimators based on As before, smoothing approaches become more interesting as the number of discretization points and the noise level increase. The empirical coverage of the smoothed estimator is lower than the linear interpolation estimator but remains slightly higher than the nominal one.
SRSWOR  Stratified sampling  

Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  97.7  10.79  10.63  10.79  10.95  97.9  6.03  5.95  6.02  6.11 
97.6  10.76  10.59  10.77  10.92  48.4  5.64  5.57  5.63  5.72  
97.6  10.76  10.59  10.77  10.92  97.6  5.89  5.82  5.89  5.97  
97.6  10.76  10.59  10.77  10.92  97.6  5.96  5.88  5.96  6.04  
97.7  10.50  10.33  10.50  10.65  97.8  5.60  5.52  5.59  5.68  
25%  lin  97.6  12.69  12.52  12.70  12.86  98.3  8.59  8.49  8.59  8.68 
97.5  12.47  12.31  12.48  12.64  58.1  6.34  6.24  6.34  6.44  
97.5  12.47  12.31  12.48  12.64  97.6  7.09  7.00  7.08  7.17  
97.6  12.47  12.31  12.48  12.64  97.8  7.10  7.01  7.10  7.19  
97.9  10.50  10.33  10.50  10.66  97.6  5.59  5.51  5.59  5.67 
SRSWOR  Stratified sampling  

Mean  1Q  Median  3Q  Mean  1Q  Median  3Q  
5%  lin  97.4  21.33  21.02  21.32  21.68  96.9  5.83  5.75  5.83  5.90 
97.4  21.29  20.94  21.30  21.60  58.1  5.69  5.61  5.69  5.77  
97.4  21.29  20.94  21.30  21.60  96.8  5.79  5.71  5.79  5.87  
97.4  21.29  20.94  21.29  21.61  96.6  5.79  5.71  5.79  5.87  
97.4  20.77  20.42  20.76  21.10  97.6  5.52  5.44  5.52  5.60  
25%  lin  98.0  13.51  13.33  13.52  13.68  95.7  7.79  7.71  7.78  7.86 
97.5  12.06  11.88  12.05  12.23  72.6  7.16  7.05  7.14  7.24  
97.5  12.06  11.88  12.05  12.23  95.0  7.53  7.46  7.53  7.60  
97.6  12.06  11.88  12.05  12.22  95.6  7.58  7.50  7.57  7.66  
97.2  10.49  10.31  10.48  10.66  97.4  5.52  5.44  5.51  5.60 
As a conclusion of this simulation study, it appears that smoothing is not a crucial aspect when the only target is the estimation of the mean, and that bandwidth values should be chosen by a cross validation criterion that takes the sampling weights into account. When the goal is also to build confidence bands, smoothing with weighted cross validation criteria lead to narrower bands compared to interpolation techniques, without deteriorating the empirical coverage. Smoothing strategies which do not take account of unequal probability sampling weights lead to empirical coverage levels that can be far below the expected ones.
5 Concluding remarks
In this paper we have used survey sampling methods to estimate a population mean temporal signal. This type of approach is extremely effective when data transmission or storage costs are important, in particular for large networks of distributed sensors. Considering noisy functional data, we have built the HorvitzThompson estimator of the population mean function based on a smooth version of the sampled curves. It has been shown that this estimator satisfies a CLT in the space of continuous functions and that its covariance can be estimated uniformly and consistently. Although our theoretical results were presented in this paper with a HorvitzThompson covariance estimator, they are very likely to hold for other popular estimators such as the SenYatesGrundy estimator. We have applied our results to the construction of confidence bands with asymptotically correct coverage. The bands are simply obtained by simulating Gaussian processes conditional on the estimated covariance. The problem of bandwidth selection, which is particularly difficult in the survey sampling context, has been addressed. We have devised a weighted crossvalidation method that aims at mimicking an oracle estimator. This method has displayed very good performances in our numerical study; however, a rigorous study of its theoretical properties remains to be done. Our numerical study has also revealed that in comparison to SRSWOR, unequal probability sampling (e.g. stratified sampling) yields far superior performances and that when the noise level in the data is moderate to high, incorporating a smoothing step in the estimation procedure enhances the accuracy in comparison to linear interpolation. Furthermore, we have seen that even when the noise level is low, smoothing can be beneficial for building confidence bands. Indeed, smoothing the data leads to estimators that have higher temporal correlation, which in turn makes the confidence bands narrower and more stable. Our method for confidence bands is simple and quick to implement. It gives satisfactory coverage (a little conservative) when the bandwidth is chosen correctly, e.g. with our weighted crossvalidation method. Such confidence bands can find a variety of applications in statistical testing. They can be used to compare mean functions in different subpopulations, or to test for a parametric shape or for periodicity, among others. Examples of applications can be found in Degras (2011).
This work also raises some questions which deserve further investigation. A straightforward extension could be to relax the normality assumption made on the measurement errors. It is possible to consider more general error distributions under additional assumptions on the moments and much longer proofs. In another direction, it would be worthwhile to see whether our methodology can be extended to build confidence bands for other functional parameters such as population quantile or covariance functions. Also, as mentioned earlier, the weighted crossvalidation proposed in this work seems a promising candidate for automatic bandwidth selection. However it is for now only based on heuristic arguments and its theoretical underpinning should be investigated.
Finally, it is well known that taking account of auxiliary information, which can be made available for all the units of the population at a low cost, can lead to substantial improvements with model assisted estimators (Särndal et al. 1992). In a functional context, an interesting strategy consists in first reducing the dimension through a functional principal components analysis shaped for the sampling framework (Cardot et. al 2010a) and then considering semi parametric models relating the principal components scores to the auxiliary variables (Cardot et al. 2010b). It is still possible to get consistent estimators of the covariance function of the limit process but further investigations are needed to prove the functional asymptotic normality and deduce that Gaussian simulationbased approaches still lead to accurate confidence bands.
Acknowledgement. Etienne Josserand thanks the Conseil Régional de Bourgogne, France for its financial support (FABER PhD grant).
Appendix
Throughout the proofs we use the letter to denote a generic constant whose value may vary from place to place. This constant does not depend on nor on the arguments . Note also that the expectation is jointly with respect to the design and the multivariate normal model.
Proof of Theorem 1.
We first decompose the difference between the estimator and its target as the sum of two stochastic components, one pertaining to the sampling variability and the other to the measurement errors, and of a deterministic bias component:
(22) 
where and are defined in (9).
Bias term.
To study the bias term in (22), it suffices to use classical results on local linear smoothing (e.g. Tsybakov (2009), Proposition 1.13)
together with the Hölder continuity (A2) of the to see that
(23) 
Hence, for the bias to be negligible in the normalized estimator, it is necessary that the bandwidth satisfy as .
Error term.
We now turn to the measurement error term in (22), which can be seen as a sequence of random functions.
We first show that this sequence goes pointwise to zero in mean square (a fortiori in probability) at a rate .
We then establish its tightness in , when premultiplied by , to prove the uniformity of the convergence over .
Writing the vector of local linear weights at point as
and using the i.i.d assumption (A4) on the we first obtain that
Then, considering the facts that by (A1), is uniformly bounded in by (A4), and exploiting a classical bound on the weights of the local linear smoother (e.g. Tsybakov (2009), Lemma 1.3), we deduce that 

(24)  
We can now prove the tightness of the sequence of processes . Let us define the associated pseudometric
We use the following maximal inequality holding for subGaussian processes (van der Vaart and Wellner (2000), Corollary 2.2.8):
(25) 
where is an arbitrary point in and the covering number is the minimal number of balls of radius needed to cover . Note the equivalence of working with packing or covering numbers in maximal inequalities, see ibid p. 98. Also note that the subGaussian nature of the smoothed error process stems from the i.i.d. multivariate normality of the random vectors and the boundedness of the for .
By the arguments used in (Proof of Theorem 1.) and an elementary bound on the increments of the weight function vector (see e.g. Lemma 1 in Degras (2011)), one obtains that