Lin Wu

Marc Bocquet

Optimal redistribution of the background ozone monitoring stations over France

Keywords: Air quality, ozone monitoring, network design

Ozone is a harmful air pollutant at ground level, and its concentrations are routinely measured with monitoring networks. The network design problem aims at determining the optimal positioning of the monitoring stations. In this study, the background stations of the French routine pollution monitoring network (BDQA) are partially redistributed over France under a set of design objectives. These background stations report ozone variations at large spatial scale comparable with that of a chemistry-transport model (CTM). The design criterion needs to be defined on a regular grid that covers France, where in general no ozone observations are available for validation. Geostatistical ozone estimation methods are used to extrapolate concentrations to these grid nodes. The geostatistical criteria are introduced to minimize the theoretical error of those geostatistical extrapolations. A physical criterion is also introduced to measure the ability of a network to represent a physical ozone field retrieved from CTM simulations using geostatistical extrapolation methods. A third type of criteria of geometrical nature, e.g. a maximal coverage of the design domain, are based uniquely on the distance between the network stations. To complete the network design methodology, a stochastic optimization method, simulated annealing, is employed in the algorithm to select optimally the stations.

Significant improvement with all the proposed criteria has been found for the optimally redistributed network against the original background BDQA network. For instance, the relative improvements in the physical criterion value range from 21% to 32% compared to randomly relocated networks. Different design criteria lead to different optimally relocated networks. The optimal networks under physical criteria are the most heterogeneously distributed. More background stations are displaced to the coast, frontiers, and large urban agglomerations, e.g. Paris and Marseilles. The ozone heterogeneous fields are not as well reconstructed from optimal networks under geostatistical or geometrical criteria as from the optimal network obtained with the physical criterion. The values of the physical criterion for the geostatistically and geometrically optimal networks show deteriorations of about 8% and 17% respectively compared to that of the physically optimal network.

Introduction

Ozone is a harmful pollutant at ground level. Its high concentrations potentially injure human health, damage vegetation and material [START_REF] Seinfeld | Ozone air quality models: A critical review[END_REF][START_REF] Pleijel | Ozone risk assessment for agricultural crops in europe: Further development of stomatal flux and flux response relationships for european wheat and potato[END_REF]. Ground monitoring networks have been deployed to evaluate the ozone concentrations. Due to its complex chemical mechanism [START_REF] Meng | Chemical coupling between atmospheric ozone and particulate matter[END_REF] and the forcing of atmospheric transport, the ozone field is heterogeneous. This heterogeneity indicates that the spatial distribution of the monitoring stations could be optimized. The optimal positioning of the ozone monitoring stations is referred to as the ozone network design problem.

Email address: Lin.Wu@cerea.enpc.fr (Lin Wu)

A typical methodology for an ozone network design problem involves a design criterion, an ozone concentration estimation method, and an algorithm for the selection of monitoring sites. The definition of the network design problem is nevertheless problem-specific. For example, a dense redundant network might be reduced to save maintenance cost [START_REF] Nychka | Design of air quality monitoring networks[END_REF][START_REF] Fuentes | Bayesian entropy for spatial sampling design of environmental data[END_REF][START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF]. Conversely a sparse network is often sought to be redistributed or augmented [START_REF] Nychka | Design of air quality monitoring networks[END_REF][START_REF] Rayner | Optimizing CO 2 observing networks in the presence of model error: results from transcom 3[END_REF].

A network station has its own spatial scale representativeness of the underlying ozone field. Background stations, which observe large spatial scale ozone concentrations comparable with ozone fields from the sim-ulations of chemistry-transport models (CTM), are employed in many applications especially for the validation of CTMs. In this paper, we will examine how a selection of the background stations in the French BDQA (Base de Données sur la Qualité de l'Air) network could be redistributed to a regular grid over France for a better performance. This is of practical concern at French operational centers 1 .

In our redistribution methodology, we use a geostatistical ozone estimation method (also called kriging) and a station selection algorithm based on simulated annealing. Both were developed and described in detail in [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF]. The main difficulty is that the design criteria have to be defined for each regular grid cell of the design domain where no observations are available for validation. The design criteria are rather subjective [START_REF] Müller | Collecting spatial data: optimum design of experiments for random fields, 3rd Edition[END_REF]Abida et al., 2008), although some may be favored by the physical context. Different criteria may lead to different optimal networks.

The simplest criterion is a measure of the distance between stations in the network, upon which geometrical criteria, e.g. the space filling design criteria, can be defined for a better geometrical coverage of the domain [START_REF] Nychka | Design of air quality monitoring networks[END_REF].

Complex criteria can be defined using geostatistical estimation methods, which are essentially based on the spatial correlation structure, to extrapolate concentrations to those regular grid cells. Geostatistical criteria, that seek to minimize the theoretical extrapolation error, make use of this correlation structure. This is closely related to the optimum experimental design theory [START_REF] Fedorov | Optimal experimental design: spatial sampling[END_REF][START_REF] Müller | Collecting spatial data: optimum design of experiments for random fields, 3rd Edition[END_REF]. Another set of popular criteria is defined with the notion of entropy, in which the reduction of uncertainty given the network observations is maximized (see [START_REF] Le | Statistical analysis of environmental space-time processes[END_REF] for a review).

In the above criteria, the dynamics of the chemistry and transport of the ozone field are not considered explicitly. In this regard, a third set of criteria, the physical criteria, can be introduced to assess the ability of a network to reconstruct a reference physical ozone field from the observations sampled at the network sites. This reference ozone field can be generated, for instance, as the simulation output of a chemistry-transport model.

It is of both theoretical and practical value to compare the optimal networks under different criteria. For instance, [START_REF] Johnson | Minimax and maximin distance designs[END_REF] demonstrate that the space filling design of minimax type is identical to the geostatistical criterion that minimizes the maximum extrapo-1 private communication from L. Rouil and B. Besagnet at INERIS lation error, if the correlations are supposed to be independent. Under Gaussian assumptions, [START_REF] Lee | On the equivalence of kriging and maximum entropy estimators[END_REF] show that the kriging and maximum entropy estimators are equivalent, which implies that in some cases the entropy criterion is identical to the geostatistical criterion.

The main contribution of this paper is to compare three types of criteria and their impact on the redistribution results, as well as to assess the performance of the optimal networks for ozone nowcasting. For instance, one could wonder whether geometrical criteria are good substitutes for physical or geostatistical criteria [START_REF] Nychka | Design of air quality monitoring networks[END_REF], and whether geostatistical criteria are good substitutes for physical criteria. Hopefully, practical instructions for the redistribution of BDQA background stations can also be learnt from the theoretical optimization results.

The paper is organized as follows. The network redistribution methodology is presented in Sec. 2, where the design criteria, the geostatistical ozone estimator, and the selection algorithm of the background stations are detailed. Section 3 presents the setup of the redistribution experiment. The redistribution results are compared and discussed in Sec. 4. Conclusions are provided in Sec. 5.

Methodology

Background stations

An ozone monitoring station provides instantaneous or averaged observations on ozone concentrations at the measuring location. According to the local scenarios, e.g. the emission rates and the meteorological conditions, the monitoring stations have their specific representativeness of a certain spatial scale. This can be roughly described by the typology of the stations. For example, in general, rural stations record ozone variations with a spatial scale larger than that of urban stations. By contrast, an industrial station may be highly influenced by the local emissions, thus its scale representativeness may be less than one kilometer. Simultaneous treatments of observations from stations of multiple spatial scales are far from straightforward [START_REF] Malherbe | Using different spatial scale measurements in a geostatistically based approach for mapping atmospheric nitrogen dioxide concentrations. Application to the French Centre region[END_REF]. Nevertheless, it is improper to strictly relate the scale representativeness to the station typology. For instance, an industrial station at rural region may behave more like a rural station, when the industrial site around is less productive or simply shut down.

A typical Eulerian chemistry-transport model computes the concentrations of a set of chemical species by solving a system of advection-diffusion-reaction equations (see [START_REF] Sportisse | A review of parameterizations for modelling dry deposition and scavenging of radionuclides[END_REF] for a review). The evolution of the species concentrations depends on many factors, e.g. the meteorological conditions, emission and deposition rates, and the chemical reactions among species. In general, the concentrations are computed at a regular grid defined in a three-dimensional model domain.

For regional applications, the model grid interval is usually from 5 to 50 km. For a successful fusion of the information from both the CTM model and observations, the monitoring stations are expected to have comparable scale representativeness. In other words, we are interested in a subset of network stations, which monitor large scale ozone field and will be referred to as background stations in this paper.

Yet there is no clear definition for the background station in the literature. We will select the background stations among the BDQA stations by their accordance with the CTM simulations over France. This accordance can be evaluated by the root mean square error and/or the correlation between the station observations and the CTM simulations. Note that the model error, which addresses the model deficiency due to the discretization (thus with unresolved smaller scales), is not explicitly taken into account.

Spatial interpolation

A geostatistical ozone estimation method is employed in our network redistribution problem. The algorithmic details can be found in [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF]. The ozone concentrations are regarded as realizations of a spatiotemporal random field Z with given means and a known spatial correlation structure.

Let G = {s 1 , . . . , s p } be the set of gauged locations of p monitoring stations, and let U = {s p+1 , s p+2 , . . . , s n } be the set of ungauged locations of given regular grid points targeted for redistribution in the French territory. The ozone concentration vector x = Z(s p+1 ), Z(s p+2 ), . . . , Z(s n ) T at ungauged sites U can be estimated by a best linear unbiased estimator (BLUE) x, which is a linear combination of the observation vector y = Z(s 1 ), Z(s 2 ), . . . , Z(s p ) T at gauged sites G:

x

= E[x] + Σ xy Σ -1 yy (y -E[y]) , (1)
where

Σ xy = E (x -E[x])(y -E y ) T , (2)
Σ yy = E (y -E[y])(y -E y ) T , (3)
are the covariance matrices. This estimator minimizes the total variance of the unbiased estimation error ǫ = xx. The covariance matrix for the estimation error is

Σ ǫǫ = E ǫǫ T = Σ xx -Σ xy Σ -1 yy Σ yx . (4)
The spatial correlations determine how the information from observation y is dispatched in the domain by the BLUE analysis. The covariance matrices Σ xx , Σ yy , Σ xy , Σ yx are computed based on an isotropic diurnal nested model:

C (h) = c 0,b + σ 2 b e -h L b , b ∈ I b , ( 5
)
where h is the distance between two sites, C (h) is the covariance of ozone concentrations between the two sites, I b = {0, . . . , n b -1} is the set of indices of bins which partition the 24 hours of a day into several intervals of equal time length, and n b is the total number of bins. Here for each bin indexed by b, L b is the correlation length, c 0,b is the background correlation for long distances, and σ 2 b is the a priori variance of the error field with subtraction of the background influence. The variance for the field is thus c 0,b + σ 2 b . For hourly bins, n b equals to 24. Note that in this formula, the ozone field is assumed to be daily stationary, but time-varying during the day. The parameters of the covariance model should be calibrated accordingly.

Three types of redistribution criteria

The network redistribution problem can be formalized by:

ξ * = argmin ξ Ψ(ξ) , ( 6
)
where ξ is a potential network configuration, and Ψ is a certain scalar criterion. Again, denote G = {s 1 , . . . , s p } as the set of locations of p monitoring stations, and U = {s p+1 , s p+2 , . . . , s n } the set of locations of the ungauged grid points targeted for redistribution. Let B be the binary set {0, 1}, ξ ∈ B n is then the vector that describes the network configuration. The i-th component ξ i of the configuration vector is 1 if site s i is included in the displaced network, otherwise ξ i = 0. Let r < p be the number of stations to be displaced, the configuration ξ must satisfy the following constraints:

n i=1 ξ i = p , n i=p+1 ξ i = r . (7)
Let V 1 = {s i l |l = 1, . . . , r; i l ∈ {1, . . . , p}} the set of station locations to be displaced, and let V 2 = {s j l |l = 1, . . . , r; j l ∈ {p + 1, . . . , n}} the set of targeted sites for redistribution, then the observation y for spatial interpolations is defined at p points from the set B ≡ (G\V 1 ) ∪ V 2 , and the estimator x is defined at npr points from the set A ≡ U\V 2 .

The network redistribution problem is closely related to the network augmentation problem, in which the network configuration ξ is defined only for U since the existing network G is fixed.

Three types of criteria can be defined for redistribution at locations where no observations are available for validation. These criteria are of geometrical, geostatistical, or physical nature. They are described in the following.

Geometrical criteria

When there is little a priori knowledge available on the statistics and the physical properties of the ozone field, space filling design methods can be used [START_REF] Nychka | Design of air quality monitoring networks[END_REF]. The criteria are uniquely based on geometrical consideration. For instance, the optimal network may be expected to uniformly cover the design domain.

For our network redistribution problem, one geometrically optimal configuration could be

ξ * = argmax ξ                    min s 1 , s 2 ∈ A s 3 , s 4 ∈ B {d(s 1 , s 2 ), d(s 3 , s 4 )}                    , (8)
where d(•, •) is the great circle distance between two points. The minimal distance of the displaced network is to be maximized, thus the resulting network would expand and cover the whole design domain.

Minimization of kriging error

The design criteria, that needs to be defined on the ungauged sites of the regular grid, can be based on the geostatistical estimation error Eq. (4) [START_REF] Cressie | Statistics for spatial data[END_REF]. In our case, the error covariance matrix Σ ǫǫ for spatial interpolation depends on the network configuration. Then one demands that the estimation error covariance Σ ǫǫ be minimal for the optimal displaced network, that is,

ξ * = argmin ξ f (Σ ǫǫ (ξ)) , (9)
where f is a function that maps a matrix into a scalar. This function f has different forms under different optimality conditions. We adopt the terminology from the optimum experimental design theory [START_REF] Silvey | Optimal design[END_REF][START_REF] Müller | Collecting spatial data: optimum design of experiments for random fields, 3rd Edition[END_REF], in which the the covariance matrix (information matrix) for the estimated parameters is referred to. Note that in our case, the covariance matrix Σ ǫǫ is defined for the estimation error as in Berliner et al. (1999).

We test the following optimality conditions:

• A-optimality: the total estimation variance is to be minimized, that is, f = Tr (Σ ǫǫ ).

• D-optimality: the confidence region of the estimation is to be minimized. This confidence region (ellipsoid under Gaussian assumptions) can be described by the principal directions of Σ ǫǫ . The volume of that ellipsoid is proportional to the determinant of Σ ǫǫ . In this case, one can choose

f = i λ i
where λ i is the i-th eigenvalue of Σ ǫǫ .

• E-optimality: similar to D-optimality but the maximal direction of the confidence ellipsoid is to be minimized, that is,

f = max λ, λ ∈ {λ i , i = 1, . . . , n -p -r}.
It is easy to adapt f to the diurnal covariance model. For example for A-optimality,

f = b Tr Σ ǫǫ,b with b ∈ I b .
Note that [START_REF] Lee | On the equivalence of kriging and maximum entropy estimators[END_REF] show that certain maximum entropy criterion, e.g. the variability-absorption design, is identical to the D-optimality criterion under Gaussian assumptions. Thus the criterion of Doptimality defined above is a representative of the class of entropy-based criteria.

The kriging methods are statistical in nature, and are not truly faithful to the ozone variability owning to chemistry and physics. The physical properties of ozone field are accounted for approximately by using diurnal covariance models. The correlations are approximated by empirical (calibrated) models. If errors in covariance parameters are important, the criteria based on kriging error (variance or covariance) is limited and may result in antithetical networks compared to those generated under the criteria aiming at estimating covariance parameters [START_REF] Zimmerman | Optimal network design for spatial prediction, covariance parameter estimation, and empirical prediction[END_REF]. However, in the context of ozone estimation, it has been found (Fig. 6 in [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF]) that the kriging is quite robust with different calibrated covariance parameters. Therefore the limitation of geostatistical criteria is less constraining for our ozone network design.

Performance of the ozone field reconstruction

Simulations of chemistry-transport models can provide forecast ozone concentration on the regular ungauged station grid [START_REF] Reynolds | Mathematical modeling of photochemical air pollution-I: Formulation of the model[END_REF][START_REF] Russell | NARSTO critical review of photochemical models and modeling[END_REF], thanks to the numerical modeling of the pollutant chemistry and physics.

It has been found in [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF] that the CTM simulations (of 0.25 • horizontal resolution) are of limited use for kriging with observations from the complete BDQA network. This is due to the mismatch of the spatial scales between the CTM and the observations. When observations are taken only from background stations, such limitation could be reduced.

Given an ozone field taken as a reference, one can reconstruct that reference ozone field from sampled ozone concentrations (of the field) at background stations using spatial interpolation methods. It is expected that a better network of redistributed background stations has better ability in representing the reference ozone field. Then, one redistribution criterion can be defined as the (hourly) root mean square error () between the reference ozone field and the reconstructed ozone field:

ξ * = argmin ξ         1 |T | × |A| k∈T ,s∈A ( xk (s) -z k (s)) 2         1 2 , ( 10
) where | • | denotes the set cardinality, k is the time index in a given period T , z k (s) is the concentration of the reference ozone field at time k and site s ∈ A, and xk (s) is the corresponding ozone estimation using the spatial interpolation formula Eq. ( 1). The values of the observation y for kriging take the reference realizations z(s) at monitoring sites in B.

This criterion is consistant with statistical indicators () used in air quality modeling and it is related to several European air quality standards.

In addition to the CTM simulations, concentration observations can be obtained by direct measuring of the ozone field. However, the utility of the observations are limited, because the measurements are in general erroneous and conducted at a given spatiotemporal scale. Data assimilation algorithms merge the information from both the CTMs and the observations aiming at a better account of the true ozone field [START_REF] Elbern | Ozone episode analysis by fourdimensional variational chemistry data assimilation[END_REF][START_REF] Wu | A comparison study of data assimilation algorithms for ozone forecasts[END_REF]. Thus, the assimilation results can also serve as a reference ozone field for reconstruction.

A third possible reference ozone field can be generated with the kriging results on the sites of U using observation from the location set G of the background stations.

In summary, we will test seven design criteria in this paper. They are of different nature and can be arranged in three sets listed in Tab. 1.

Optimization using simulated annealing

In the network design context, the combinatorial design problem Eqs. (6, 7) can be solved by the simulated annealing algorithm [START_REF] Van Groenigen | Constrained optimization of spatial sampling using continuous simulated annealing[END_REF]Abida et al., 2008;Abida and Bocquet, 2009;[START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF]. In a classical annealing, the optimization process can escape from a local minimum ξ (i) to a new configuration ξ (i+1) with an acceptance probability [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]:

P(ξ (i) , ξ (i+1) , τ) = min             1, e - Ψ(ξ (i+1) ) -Ψ(ξ (i) ) τ             , ( 11
)
where τ is a global parameter which is an analog of temperature. Often τ is initially high, and the iterative process probes large-scale variation of Ψ. When τ decreases according to certain cooling schedule, the iterations search for finer variations. By carefully choosing the cooling schedule, the global minimum can be approached to some precision which can be arbitrarily small. We employ the geometric cooling schedule,

τ (k+1) = ατ (k) , (12)
where α (0 < α < 1) is a decreasing factor, and usually k coincides with i.

The new candidates are chosen from the neighborhood of ξ (i) via a flipping procedure. The configuration ξ (i) is divided into two parts: the first p components for the monitoring stations, and the last np components for the targeted points. For each part, we randomly flip one component's value from one to zero, and randomly flip another component's value from zero to one. By this way, we firstly randomly choose one different station to be displaced, then randomly move this station to a different targeted point. Note that the constraint Eq. ( 7) is automatically satisfied. In practice, the tuning of the parameters values, especially for the initial and final temperatures, is necessary to obtain a satisfactory solution.

Experiment setup

CTM simulation

For this study, the P/P3D model [START_REF] Boutahar | Development and validation of a fully modular platform for numerical modelling of air pollution: POLAIR[END_REF][START_REF] Mallet | Technical Note: The air quality modeling system Polyphemus[END_REF][START_REF] Sartelet | Simulation of aerosols and gas-phase species over Europe with the Polyphemus system: Part I-Model-to-data comparison for 2001[END_REF] [START_REF] Troen | A simple model of the atmospheric boundary layer; sensitivity to surface evaporation[END_REF]) (in the unstable boundary layer); 9. boundary conditions: gas and aerosol concentrations from the general circulation model LMDz [START_REF] Hourdin | The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection[END_REF]; 10. numerical schemes: a first-order operator splitting, the sequence being advection-diffusionchemistry; a direct space-time third-order advection scheme with a Koren-Sweby flux-limiter; a second-order order Rosenbröck method for diffusion and chemistry [START_REF] Verwer | Numerical time integration for air pollution models[END_REF]. Ozone concentrations peak during summer, which is the most risky scenario for human health and crop production. The CTM simulations covers summer 2005 from 20 May at 0000 UTC to 1 September at 0000 UTC. The initial conditions are interpolated based on the coarser LMDz simulations. The first 12 days are the spin-up period, and the simulations after spin-up are used in the network redistribution problem.

Background observations

The BDQA (Base de Données sur la Qualité de l'Air) network is dense for regional applications (information available at http://www.atmonet.org). There are 678 stations within France. These stations (see Fig. 1) are located in typologically different areas, such as urban districts and regional areas of cities, industry sites, and heavy traffic roads. The background stations are observers for the large scale chemistry-transport phenomena comparable with the CTM simulations (see Sec. 2.1). For each BDQA station, we collect its observations of hourly ozone mean concentrations from 1 June at 0000 UTC to 1 September at 0000 UTC. Then the correlations and  between those observations and the P-/P3D simulations for this given period are computed. In this paper, background stations are defined as those with correlation larger than 0.70 and with  smaller than 25 µg m -3 . According to this definition, for summer 2005, there are 190 background stations (shown in Fig 1). Compared with the full BDQA network, most of the rural stations are picked up as background stations. The majority of the retained sites are rural and urban stations. There are also five industrial and one traffic stations in the set of background stations. The stations in south-east and the France-Spain border tend to be excluded in the network of background stations, because at these mountainous areas (Alpine territory), it is difficult to represent a large scale ozone field. These statistics are consistent with the choice of BDQA background stations from INERIS5 .

Assimilation using optimal interpolation

The information from the CTM simulations and BDQA observations can be combined by data assimilation algorithms to generate a more realistic ozone field.

At one time step, the ozone concentrations computed by the CTM simulations (denoted as c b ) are adjusted by the observation vector o, so that the error variance of the ozone estimation c is minimized. Under the Gaussian assumptions, a BLUE formula for optimal interpolation (OI) reads:

c = c b + BH T (HBH T + R) -1 (o -H(c b )) , ( 13
)
where H is the observation operator that maps the ozone concentration vector to the observation vector, H is the linear form of H, B is the covariance matrix of the errors in c b , and R is the (assumed diagonal) covariance matrix of the errors in o. The matrix B is parameterized in Balgovind form, that is, the error correlations are isotropic and described by the Balgovind function:

f (h) = (1 + h L B )e -h L B σ 2 B ,
where h is the distance between two locations, L B is a characteristic length and σ 2 B a background error variance. In this study, the parameters L B and σ 2 B are set to a priori values as 1 • and 400 (µg m -3 ) 2 respectively. The observational error variance takes value of 100 (µg m -3 ) 2 .

It has been shown in [START_REF] Wu | A comparison study of data assimilation algorithms for ozone forecasts[END_REF] that, although relatively simple compared with other advanced assimilation algorithms, OI is an effective assimilation method for ozone estimation. In this study, only isotropic correlations are considered. Anisotropic considerations could bring a positive impact, but the improvements are not significant (Fig. 7 in [START_REF] Blond | Three-dimensional ozone data analysis with an air quality model over the Paris area[END_REF]). Furthermore, only observations from background stations are used in this paper, which lessens the need to introduce an anisotropic correlation model.

Kriging over the full domain

A reference ozone field can also be generated by kriging over the full domain U using Eq. (1). We use a constant mean. The covariance model in use is with one hour bins and calibrated to the BDQA background station observations (detailed in Sec. 4.1). The period for kriging is also taken from 1 June to 1 September in 2005.

Results and discussions

Covariance models and kriging performance

For a given time window T , the covariance for the ozone concentrations at two sites s i , s j within the time period of a given bin can be estimated by

C(s i , s j ) = 1 N t N t t=1 Z t (s i ) -Z(s i ) Z t (s j ) -Z(s j ) , (14) with
Z(s i ) = 1 N i N i t=1 Z t (s i ) , (15)
where N i is the number of observations at site s i , N t is the number of mutually available observation pairs for site pair (s i , s j ), and Z t (s i ) denotes the random variable for the ozone field at time index t and site s i .

When s i , s j run on all the available sites, a cloud of covariance values against the distance can be obtained. The covariance cloud can then be averaged within continuous regions

T (h i ) = [h i -L T /2, h i + L T /2], for h i = i × L T -L T /2, i ∈ N.
Here L T is set to 30 km in this paper for the regional application. By this way, the curves of the regionalized covariances can be plotted. We omit the figures about the covariance clouds and the regionalized covariance curves, because they are very similar to Fig. 2 and Fig. 3 in [START_REF] Wu | A comparison study of data assimilation algorithms for ozone forecasts[END_REF].

Let Θ = [c 0,b , σ 2 b , L b ] be the vector of unknown parameters for the nested covariance function C(•) in Eq. ( 5). The parameter Θ is determined by solving the ordinary least-square fitting problem

Θ * = argmin Θ N h i=1 C(h i ) -C(h i ) 2 (16)
where N h is the total number of the tolerance regions, h i is the center of the i-th region, and C(h i ) is the corresponding regionalized covariance. For each bin, the parameter vector Θ is calibrated using the Levenberg-Marquardt algorithm. For observations at available background BDQA stations during summer 2005, the fitting results for one-hour bins are shown in Fig. 2. The correlation length L and the variance for the field c 0 + σ 2 are smaller at night and peak at noon. This is also similar to the fitting results with observations from all BDQA stations but during a shorter period of one month [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF].

When Z t (s i ) is sampled from the CTM simulations or the data assimilation results, covariance clouds and regionalized covariances can be obtained in a similar way. Accordingly the covariance models can be calibrated for hourly bins. The fitted parameters for these hourly covariance models are shown in Fig. 3. These parameters clearly demonstrate a diurnal cycle. The parameters of the covariance model for analyzed ozone field is a compromise between the covariance models for BDQA observations and CTM simulations. Its correlation length shows a similar pattern to that of the covariance model fitted to BDQA observations, and its variance is a fusion of the two other covariance models. We evaluate the reconstruction  Eq. ( 10) with reference ozone field set to the CTM simulation (criterion PHY-S in Tab. 1) for a large set of randomly displaced networks. These networks are generated by randomly displacing a given number of background stations to the regular ungauged grid points within French territory. The reconstruction  are shown in Fig. 4. The worst reconstruction  is less than 5 µg m -3 . This is far inferior to the CTM model error for hourly ozone concentrations (typically more than 20 µg m -3 ). This indicates that the kriging method is satisfactory in reconstructing the ozone field generated by CTM simulations. The covariance model calibrated to BDQA observations is employed for the above evaluation. Two other CTM related covariance models (see Fig. 3) are also tested for the same set of random networks. It is found that the reconstruction  is not sensitive to the covariance models. The maximal relative  difference is 2% for these random networks. For the background BDQA network (see Fig. 1), the reconstruction s for the three covariance models are 4.61 , 4.57 and 4.57 µg m -3 .

Results of a reference redistribution optimization

Before carrying out sensitivity studies, we first define a reference redistribution optimization. In this reference algorithmic setting, the parameters for the hourly covariance model are fitted to the observations from BDQA background stations. There are only 40 stations to be displaced to the regular grid of ungauged sites over France. A geometric schedule is employed with annealing rate set to 0.9999. The criterion for redistribution is chosen to be the reconstruction  of the reference ozone field generated by CTM simulations (PHY-S in Tab. 1).

The optimal reconstruction  for the reference setting is 3.121 µg m -3 . A relative 28% improvement in the reconstruction  is obtained through optimal redistribution compared with the case of the background network without redistribution. The resulting map of the redistributed network is shown in Fig. 5a. The clustered stations tend to be moved to the regions with sparse network coverage, e.g. the south-east France and the frontier area between France and Spain. In these areas, there are lacks of background stations. Target locations in these large regions are uniformly distributed.

Sensitivity of the design to algorithmic settings 4.3.1. Discrepancy between two networks

The discrepancy between two different networks needs a quantitative definition, so that the interpretation of the optimization results could be better approached. A geometrical distance between two networks, say B a and B b , has been introduced by [START_REF] Nychka | Design of air quality monitoring networks[END_REF] as Another global measure on the discrepancy between B a and B b for a given region D is to compute the difference of their coverage for each site s ∈ D [START_REF] Saunier | Model reduction via principal component truncation for the optimal design of atmospheric monitoring networks[END_REF]. The coverage of a network B at s can be estimated by the number of the network stations located within a circle centered at s of radius R D . Let us denote this coverage number as λ B (s). The coverage difference

∆(B a , B b ) = max
δλ B a ,B b (s) = λ B a (s) -λ B b (s), s ∈ D ( 18
)
provides a smooth discrepancy between two networks B a and B b . The map of this smooth discrepancy describes the spatial difference (local augmentation or depletion of stations) between networks. In our case, D is the regular ungauged grid over France, and the probe radius R D is set to 75 km. Eventually an average of the absolute coverage difference

Λ(B a , B b ) = 1 |D| s∈D |δλ B a ,B b (s)| (19)
serves as another scalar measure of the difference between two networks.

Sensitivity to simulated annealing and covariance models

In this section, we test the sensitivity of the redistribution results to the algorithmic settings. Five additional redistribution experiments are performed; each one with a single alternating parameter different from the reference algorithmic setting in Sec. 4.2. These alternative parameters are either about simulated annealing or about the nested covariance model.

Figure 5 describes the setting details (read the caption) and shows the resulting redistributions with both reference and the five alternative algorithmic settings. The optimal reconstruction  for the five alternative algorithmic settings (b-f) are 3.139 , 3.159 , 3.126 , 3.152 and 3.132 µg m -3 respectively. The maximal relative difference in reconstruction  for these alternative algorithmic settings against the reference setting (3.121 µg m -3 ) is 1.2%. The maps of the redistributed networks with all the six algorithmic settings show a very similar pattern. The resulting redistribution is only slightly sensitive to the covariance model and the annealing setting. Some target locations for the displacement are even identical near country borders.

In the following, if not mentioned, the optimizations are performed using the covariance model fitted to BDQA data (Fig. 2) with the slow annealing rate (0.9999). 3a), and in the setting (f) the covariance model is fitted to the assimilation results (Fig. 3b). The circles show the remaining stations, and the cross signs show the locations of the stations chosen to be reallocated. The squares indicate the new locations of these chosen stations.

Sensitivity to the redistribution criteria

The maps of the redistributed networks under the seven criteria listed in Tab. 1 are shown in Fig. 6. The optimal network under the geometrical criterion MAXMIN is uniformly distributed to cover the French continental territory. The optimal networks under the geostatistical criteria are less uniform because of smaller correlations in between distant stations. The optimal networks generated under criteria STAT-A and STAT-E are very similar. By contrast, the optimal network under criterion STAT-D are more evenly distributed than its two counterparts, since the volume of the confidence ellipsoid defined by the estimation error covariance matrix Σ ǫǫ are minimized, whereas the other two criteria assign more weights to the main direction of the confidence ellipsoid.

The optimal networks under the physical criteria are more heterogeneous than the above two catalogues of optimal networks. It is reasonable because the reference ozone field itself is heterogeneous. More stations are allocated around coast, frontiers, and large urban agglomerations, e.g. Paris and Marseilles, so that the high and uncertain ozone concentrations at these regions can be better represented. The assimilation of concentration observations would render the reference ozone field more realistic and increase the ozone variation around the uncertain regions compared with the plain CTM simulations. As a result, the clusterings of stations at the uncertain regions are more evident in the optimal network with assimilation results as reference field (criterion PHY-A; Fig. 6a) than that with CTM simulations as reference field (criterion PHY-S; Fig. 5a).

When the reference ozone field is generated via kriging using ozone observations from all the background stations (criterion PHY-K; Fig. 5b), the new locations remain close to those of the original background stations, which are the best candidates in reconstructing this reference field.

We list the performance of the resulting optimal networks under the seven criteria in Tab. 2. For any criterion, the network with the best performance is the optimal network obtained under that criterion. The optimal network under criterion PHY-K remains close to the original background network and performs the worst under other criteria. The geometrically optimal network under criterion MAXMIN is the most uniformly distributed and performs second worst under other criteria.

The optimal networks for the other two physical criteria PHY-S and PHY-A perform similarly under other criteria. The optimal networks for the three geostatistical criteria also perform similarly. The optimal networks under geostatistical or geometrical criteria perform poorer for the representation of the heterogeneous ozone field. The values of the reference physical criterion PHY-S for the geostatistically and geometrically optimal networks increase by about 8% and 17% respectively with respect to that of the physically optimal network. Therefore, for complicated design objectives, e.g. those related to the heterogeneity of the ozone field, the physical criterion may be more appropriate than the geostatistical or geometrical criteria.

The geometrical distance Eq. ( 17) and the average coverage difference Eq. ( 19) between these seven optimal networks are computed and listed in Tab. 3. The optimal networks under criteria MAXMIN and PHY-K have larger distances (or differences) with other optimal networks. This is reasonable since the two optimal networks are either the most uniform or the most heterogeneous among the seven optimal networks.

Figure 7 shows the map of the smooth discrepancy Eq. ( 18) between the optimal networks under criteria PHY-S and MAXMIN. It can be clearly observed that more stations are displaced to the south-east coast and the frontier area near Spain for the optimal network under PHY-S. Note that there are slightly more stations around Paris for the optimal network under criterion MAXMIN. The number of stations to be displaced for this case is set to 40, which is not large enough to break the station clustering around Paris.

Sensitivity to the number of stations to be displaced

In this section, we evaluate the impact of the number of stations to be displaced on the design results. The cri- terion is chosen to be PHY-S. The reconstruction s of the optimal networks against the number of displaced stations are shown in Fig. 4. The relative improvements in reconstruction  range from 21% to 32% as compared to the random networks. As expected, with more stations displaced, one observes improvement in the reconstruction . However, there are no obvious improvement beyond 100 displaced stations. On the contrary the improvement seems to diminish.

This might be explained by the differences in the sites' locations between the BDQA network from which the stations are drawn to be moved, and the regular grid of targeted sites for relocation. Since there is a minimal distance between these grid points, there will also be a minimal kriging error when interpolating from gauged grid points to ungauged points on the grid. The error could be smaller than this minimal threshold if gauged sites are set within grid cells. In fact, a significant fraction of the BDQA stations are located in the grid cells, and a few of them are crucial for the interpolation and may have locations which are better candidates than the targeted grid points for kriging. Beyond 100 relocated stations, it seems that part of these crucial stations are being relocated to worse sites (grid points), which starts degrading the global performance.

The maps of optimally displaced networks with different number of displaced stations are shown in Fig. 8. When a few stations are displaced, the new locations tend to be in the regions of heterogeneous and uncertain ozone concentrations, e.g. the coast, the frontiers, and the large urban agglomerations. When more stations are displaced, the new locations result from a balance between the ozone heterogeneity and the coverage of the design domain.

Population density as an additional constraint

Other spatial constraints may also be taken into account in addition to the three sets of criteria introduced in this study, such as the local population density and the cost of the stations. This can be considered as a multiobjective network design, but detailed investigations on this matter are beyond the scope of this paper. However, we conduct a preliminary study by incorporating the local population density in our criteria.

The physical criterion Eq. ( 10) is adjusted to take into account the local population density at any point in the domain U: where w(s) is the population density at grid point s, and W is the total population in the complete design domain U or in the estimation domain A ⊂ U. With this new criterion, the regions with large population are given more importance in the optimization. This addresses the fact that the emissions are closely related to the human activities and that the exposure to the air pollution at populated regions would have more social and economic impact. Whether the population is summed over the complete or the estimation domain makes little difference since the size of the set of displaced stations V 1 is much smaller than the size of the set of estimation points A.

ξ * = argmin ξ         1 |T | × |A| × W k∈T ,s∈A w(s) (x k (s) -z k (s)) 2         1 2 , (20)
The map of the optimal redistributed networks under this new criterion is shown in Fig. 9, in which the reference ozone field is taken to be the CTM simulation. It is clear that, under this criterion, the stations tend to be displaced to the regions with large population ( Île-de-France and Vallée-du-Rhône regions). Note that there are still many stations moved to the regions with lower population density. This is necessary for a reasonable reconstruction of the reference ozone field.

Conclusion

We have developed a methodology for the relocation of the BDQA background stations that monitor large spatial scale ozone properties comparable with the scales probed by an Eulerian chemical and transport model. This methodology relies on kriging methods for ozone estimations, on a set of redistribution criteria, and on a selection algorithm based on simulated annealing. The kriging performs very well in the reconstruction of the reference field generated by CTM simulations.

The monitoring network redistribution problem turns out to be quite different from the reduction problem studied by [START_REF] Wu | Optimal reduction of the ozone monitoring network over France[END_REF] because it requires estimating pollutant concentrations on sites where no observations are available for validation.

Three sets of design criteria have been defined in the design domain (a regular grid) where no observations are available for validation. These criteria are of geostatistical, physical or geometrical nature. Significant improvements on the performance for all those criteria have been achieved through optimal redistribution of the original BDQA background stations. For instance, the relative improvements in the value of a reference physical criterion range from 21% to 32% compared to the randomly relocated networks. Moreover, we have shown that this improvement saturates beyond 100 relocated stations. The optimal redistribution results have been found to be very little sensitive to the algorithmic setting, e.g. the annealing rate and the parameters of the covariance model.

The optimal redistributed networks under different criteria have been shown to be quite different. For example, the optimal networks for the physical criteria are more heterogeneous (adapting to the heterogeneous ozone reference field) than those for the other two sets of criteria. The values of the reference physical criterion (root mean square error of the reconstructed hourly concentrations) increase by about 8% for the geostatistically optimal networks and by 17% for the geometrically optimal networks, as compared to the minimum of this reference physical criterion obtained for the physically optimal network. Therefore, for complicated design objectives, e.g. those related to the heterogeneity of the ozone field, the physical criterion may be more appropriate than the geostatistical or geometrical criteria. Note that the design efficiency (computational aspect) is not considered as a factor for choosing design criteria in this paper.

For all the optimally redistributed networks, stations from the clusters of background stations tend to be moved to the regions with sparse network coverage. For the optimal networks under physical criteria, more background stations are displaced to the coast, frontiers, and large urban agglomerations, e.g. Paris and Marseilles. The optimal network under the geometrical cri-terion is the most uniformly distributed. By contrast, the optimal networks under geostatistical criteria are less uniform because of small correlations in between distant points.

In this study, errors of the model and of the observation were not taken into account. However, by selecting the background stations, we have reduced the representativity error, which is part of model error. These error terms are absent in our design criteria. An ongoing subject is to incorporate the posterior error covariance matrix, which is a natural output of the data assimilation scheme [START_REF] Wu | A comparison study of data assimilation algorithms for ozone forecasts[END_REF], into the design criteria.

The topographic constraints in complex mountain/littoral areas make them difficult to be associated with large scale ozone patterns. This renders the interpretation of the redistribution of stations to these areas less transparent. The optimal target locations in these areas could be better evaluated with nested higher resolution CTM models. The use of observations in these complex area will be a key point to better represent more complex ozone field at finer spatial scales.

We have also attempted to extend this study to other pollutants, such as PM 10 , or NO 2 . However the absence of clear long-range correlations, or that are at best not as prominent as for O 3 , gives little hope of a straightforward generalization of the schemes presented here that rely on geostatistical interpolation.

As a step forward, the augmentation of a sparse network is an affordable issue. The methodology of network redistribution can be adapted to this context, since both need the criteria defined at locations where no observations are available. The network design based on multiple criteria is another topic worth further investigations [START_REF] Reed | Simplifying multiobjective optimization: An automated design methodology for the nondominated sorted genetic algorithm-II[END_REF]. The network design in this case will result from a balance of many criteria, e.g. the proposed criteria in this paper, the constraint of population density, and the minimization of the cost of the network construction or maintenance.

Acknowledgements. This research was supported by the French Research Network on Sustainable Development (Réseau R2DS Île- de-France 2007-21)

  The model domain covers France (metropolitan area, [5.25 • W, 41.75 • N] × [9.25 • E, 52.25 • N]) with a 58 × 43 grid of 0.25 • horizontal resolution (see Fig. 1). The altitude is divided into 9 vertical layers. The tops of the vertical layers are 30 m, 150 m, 350 m, 630 m, 975 m, 1360 m, 1800 m, 2270 m and 2780 m respectively. The top layer is high enough to enclose the planetary boundary layer. A time step of 300 s is used.

Figure 1 :

 1 Figure 1: Map of the background BDQA Network. The circles indicate the locations of the background stations, and the cross signs show the locations of the BDQA stations that are not retained as background stations. The regular tiny points are the targeted candidate locations for redistribution. These points are the P/P3D model grid-cell centers within the French territory.

Figure 2 :

 2 Figure2: The calibrated parameter values with respect to bin indices for one-hour bins, using the ozone concentrations from the background stations. The unit for the background correlation c 0 and the a priori variance σ 2 is in (µg m -3 ) 2 , and the unit for the correlation length L is in km.

Figure 3 :

 3 Figure 3: The calibrated parameter values with respect to bin indices for one-hour bins. The unit for c 0 , σ 2 is in (µg m -3 ) 2 , and the unit for L is in km. The upper panel (a) shows the results using the reference simulation field. The lower panel (b) shows the results using the reference assimilation field.

Figure 4 :

 4 Figure4: The reconstruction s for the kriging with randomly and optimally displaced networks. The reference ozone field is generated by CTM simulations. For each given number of stations to be displaced, 100 random networks are generated for evaluation. The error bar shows the standard deviation in reconstruction  for these random networks. The centers of the error bar are the means of the reconstruction s. The squares shows the reconstruction s of optimal networks against the given number of stations to be displaced. The algorithmic setting for optimization is detailed in Sec. 4.2.

  s a ∈B a min s b ∈B b d (s a , s b ) . (17) where d(•, •) is the great circle distance. The distance ∆(B a , B b ) is a scalar that roughly measures a geometrical difference between two networks B a and B b .

Figure 5 :

 5 Figure 5: The maps of the redistributed networks obtained with the reference and alternative algorithmic settings, in which 40 background stations are displaced. Each alternative setting (b-f) has one changed factor compared with the reference setting (Sec. 4.2). The setting (b) employs a faster annealing rate 0.9993. The setting (c) chooses a different seed for the random number generation needed in simulated annealing. The setting (d) has a different initial network from which start the iterations. The parameters of the covariance model in the setting (e) are set to those fitted with the CTM simulations (Fig.3a), and in the setting (f) the covariance model is fitted to the assimilation results (Fig.3b). The circles show the remaining stations, and the cross signs show the locations of the stations chosen to be reallocated. The squares indicate the new locations of these chosen stations.

Figure 7 :

 7 Figure 7: Map of smooth discrepancy between optimally redistributed networks under criteria PHY-S and MAXMIN.

Figure 6 :

 6 Figure 6: The maps of the redistributed networks under different criteria. Here 40 background stations are displaced. The settings for the covariance model and simulated annealing are same as those in the reference algorithmic setting (Sec. 4.2). The criteria (a-f) are detailed in Tab. 1 in Sec. 2.3. The circles show the remaining stations, and the cross signs show the locations of the stations chosen to be reallocated. The squares indicate the new locations of these chosen stations.

Figure 9 :

 9 Figure 9: The maps of the optimally redistributed network under the reference criterion (reconstruction performance of the CTM simulations) adjusted by the population density. The population density is interpolated to the grid of the design domain. The dark regions have larger populations. The maximal population density is 1.93 million per grid cell at Paris. The logarithm of the population density data are plotted for a clear contour shape. Here 40 background stations are displaced. The circles and the cross signs show the remaining and displaced stations respectively. The squares indicate the new locations.

Figure 8 :

 8 Figure 8: Maps of optimally displaced networks. The number of displaced stations are 20, 60, 100, 140 respectively. The circles show the remaining stations, and the cross signs show the locations of the stations chosen to be reallocated. The squares indicate the new locations of these chosen stations.

  . The authors acknowledge the help from O. Saunier and O. Isnard from IRSN with the population density data, and useful discussions with scientists at INERIS and AIR-PARIF. Abida, R., Bocquet, M., 2009. Targeting of observations for accidental atmospheric release monitoring. Atmos. Env. 43, 6312-6327. Abida, R., Bocquet, M., Vercauteren, N., Isnard, O., 2008. Design of a monitoring network over france in case of a radiological accidental release. Atmos. Env. 42, 5205-5219. Berliner, L. M., Lu, Z.-Q., Snyder, C., August 1999. Statistical design for adaptive weather observations. J. Atmos. Sci. 56, 2536-2552.

PSU/NCAR mesoscale model

Global Land Cover Facility

Co-operative Programme for Monitoring and Evaluation of the Long-range Transmission of Air Pollutants in Europe

personal communication from L. Malherbe


Roman Gaignaire

Riccardo Scorretti

Ruth V Sabariego

Christophe Geuzaine
email: cgeuzaine@ulg.ac.be

Stochastic Uncertainty Quantification of Eddy Currents in the Human Body by Polynomial Chaos Decomposition

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The finite element method (FEM) can be used for computing the electric field induced into the human body by extremely low frequencies (ELF) fields [START_REF] Scorretti | Nicolas Numerical dosimetry of currents induced in the human body by ELF magnetic fields[END_REF]. To this aim, the electrical properties of the different tissues (conductivity, permittivity) are required. Unfortunately, these properties are subjected to uncertainty, and their determination is still a matter of discussion. Moreover, most of available data has been obtained from measurements on different animals and extrapolated to humans. Therefore, it is interesting to model and quantify the effect of these uncertainties on the electromagnetic fields induced into the human body. A good statistical characterisation could be obtained by applying a Monte-Carlo (MC) method. It is very time-consuming though (some days in our model). In this work, we use non-intrusive probabilistic algorithms, namely a chaos polynomial (CP) approach. Assuming that only the variances of conductivities are finite [START_REF] Gaignaire | Stochastic Uncertainty Quantification of the Conductivity in EEG Source Analysis by Using Polynomial Chaos Decomposition[END_REF], [START_REF] Gaignaire | Current calculation in electrokinetics using a spectral stochastic finite element method[END_REF], this kind of method allows to obtain a complete characterization of the induced field in a probabilistic dimension with a reduced computational cost, a few hours in comparison to MC simulations.

II. INCORPORATION OF STOCHASTIC UNCERTAINTY

A. Deterministic framework

In the deterministic context of ELF magnetic fields, it is generally assumed that the reaction field due to eddy currents in living tissues is negligible. This allows to develop specific formulations, where the computational domain Ω is limited to the human body. In this work we use the φa FE formulation with the electric scalar potential φ to determine and the magnetic vector potential a (magnetic flux density b = curl a) known a priori, as the reaction field is negligible.

This work has been supported in part by the Belgian Science Policy (IAP P6/21) and the Belgian French Community (ARC 09/14-02).

By weakly imposing div j = 0 with the boundary condition n • j| ∂Ω = 0, the φa formulation reads [START_REF] Scorretti | Nicolas Numerical dosimetry of currents induced in the human body by ELF magnetic fields[END_REF]:

(σ ∂ t a, grad φ ′ ) Ω +(grad φ, grad φ ′ ) Ω = 0 , ∀φ ′ ∈ H(grad ,Ω) .
(1) As human model, we consider the the phantom ZOL [START_REF] Scorretti | Nicolas Numerical dosimetry of currents induced in the human body by ELF magnetic fields[END_REF] available from the Visible Human project [START_REF] Ackerman | The visible human project[END_REF]. The electrical conductivity values σ are taken from [5]. The flux density source is a 50 Hz uniform vertical field b z = 1 mT. In the literature concerning the protection against ELF fields, three scalar values are generally computed from the current density j for each organ: the (spatial) average j A , the maximum value j M , and the 99% percentile [START_REF] Dimbylow | Induced current densities from low-frequency magnetic fields in a 2 mm resolution, anatomically realistic model of the body[END_REF]. In particular, the ICNIRP recommendations focus on j M . Similar definitions exist for the electric field e.

B. Uncertainties

In this paper, we are mainly concerned with the fields induced in the head. The conductivities in the brain σ B (ω), and of the cerebellum σ C (ω) are modeled within a probabilistic framework, by assuming that they are random variables. Therefore, the average j A , the maximum j M or the 99% percentiles j P of the induced fields are random as well. In particular, by applying the maximum entropy principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] we model (arbitrarily) σ B (ω) and σ C (ω) (in S/m) as independent random variables, uniformly distributed between the 1/3 and three times the values in [5], i.e.

σ B (ω) ∼ U ([0.0178 ; 0.160]) , (2)
σ c (ω) ∼ U ([0.0317 ; 0.286]) . (3)

C. The Non -Intrusive Approach

As the conductivities of the brain and the cerebellum are two independent random variables of finite variance, we can expand them as a truncated series of order p in in the twodimensional Hermite polynomials of a random gaussian vector ξ(ω) = (ξ 1 (ω), ξ 2 (ω)), known as Hermite chaos polynomials [START_REF] Gaignaire | Current calculation in electrokinetics using a spectral stochastic finite element method[END_REF]:

σ B (ω) ≈ Pin i=0 σ Bi ψ i (ξ(ω)) , (4)
σ C (ω) ≈ Pin i=0 σ Ci ψ i (ξ(ω)) , (5)
where σ Bi and σ Ci are scalar values that depend on the probabilistic law of the conductivities, P in = C pin 2+pin is the number of two-dimensional polynomial of order smaller than p in , and ψ i is the i th two-dimensional Hermite polynomial.

The stochastic problem is solved by means of a CP decomposition of both the conductivity and the induced fields [START_REF] Gaignaire | Current calculation in electrokinetics using a spectral stochastic finite element method[END_REF]. There is only an assumption on the conductivities, with no a priori hypothesis on the probabilistic distribution shape of the average j A , the maximum j M or the 99% percentiles j P of the induced fields.

Those induced fields -the average current density in the brain j A (ω) = j A (ξ(ω)) -are calculated by the FEM from any pair of values (σ B (ξ(ω)), σ C (ξ(ω))). The average density belongs to a space that can be spanned by the polynomials ψ(ξ(ω)) and can thus be written as a truncated series to an order p out :

j A (ω) = Pout i=0 j A m ψ m (ξ(ω)) . ( 6
)
The unknown real coefficients j A m are computed taking into account the orthogonality properties of the Hermite polynomials:

j A m i = E[j A (ω)ψ m (ξ(ω))] E[ψ m (ξ(ω)) 2 ] , (7)
with E[•] the mathematical expectation. The denominator can be computed analytically. A Hermite Gauss integration scheme with d points is used to determine the numerator, which is an integral [START_REF] Gaignaire | Current calculation in electrokinetics using a spectral stochastic finite element method[END_REF]:

E[j A (ω)ψ m (ξ(ω))] ≈ d i=1 ... d j=1 w i,j (j A ((t 1 , t 2 ) i,j ))ψ m ((t 1 , t 2 ) i,j ), (8)
where (t 1 , t 2 ) i,j is the i, j-th Gauss point, and w i,j the associated weight in the two-dimensional Cartesian rule. Therefore, the deterministic problem must be solved d 2 times, with the conductivity evaluated through (4) and (ξ 1 (ω), ξ 2 (ω)) = (t 1 , = t 2 ) i,j , i, j = 1, . . . , d.

III. RESULTS AND DISCUSSION

The non-intrusive polynomial chaos decomposition methods are governed by three parameters: p in determines the accuracy of the approximation made on the input random variables; p out is the order of truncation of the induced fields (maximum, average and percentiles of the electric field or current density) and d is the number of integration points. We adopt p in = 16, p out = 6 or 10, d = 14.

The probabilistic density of the average j A , the 99% percentile j P and the maximum j M current density in the brain are depicted in Fig. 1 for p out = 6 and 10. The convergence of the non-intrusive method can be clearly observed. The mean and standard deviation of j A are 1.2 10 -4 A/m 2 and 4.3 10 -5 A/m 2 , respectively. The probabilistic density of j A presents the sharpest peak among the three considered quantities. Note that the values higher than the mean are more probable than the rest. With regard to the 99% percentile j P in the white and grey matter of the brain, the mean is 1.9 10 -4 A/m 2 with a standard deviation of 6.8 10 -5 A/m 2 . It is more dispersed around the mean than the j A . Once again the right tail is slightly more probable than the left one. Finally, the maximum j M in the brain presents a mean of 3.2 10 -4 A/m 2 and a standard deviation of 1.3 10 -4 A/m 2 . This last parameter is the most spread (the standard deviation and the mean have the same order of magnitude and the probabilistic density is the largest). It is clear that the probability to be above 3 10 -4 A/m 2 is nearly 0 for j A and j P . It is not negligible anymore for j M . Note that the maximum value of computed fields may be (strongly) influenced by the FE mesh quality. Further details on the non-intrusive probabilistic method will be provided in the full paper. Special attention will be paid to the influence of p in , d and p out and some probability threshold. Results on the induced electric field and in the cerebellum will be presented as well.

Fig. 1 .

 1 Fig.1. Probabilistic density of the averaged (j A ), 99-percentile (j P ) and maximum (j M ) current density in the brain


Christopher S Withers

Saralees Nadarajah

Accurate tests and intervals based on non-linear cusum statistics

Keywords: Approximations, Cusums, Edgeworth-Cornish-Fisher, Non-linear statistics

Suppose that we are monitoring incoming observations for a change in mean via a cusum test statistic. The usual nonparametric methods give first and second order approximations for the one-and two-sided cases. Withers and Nadarajah [Statistics and Probability Letters, 79, 2009, 689-697] showed how to improve the order of these approximations for linear statistics. Here, an extension is provided for non-linear statistics with non-zero asymptotic variances. This involves development of two calculi analogous to that of von Mises.

Introduction and Summary

Suppose that we are monitoring incoming observations for a change in mean via a cusum test statistic. The usual nonparametric asymptotic methods give a first order approximation in the one-sided case and a second order approximation in the two-sided case, where by Ith order we mean an error of magnitude n -I/2 for n the sample size. [START_REF] Withers | Accurate tests and intervals based on linear cusum statistics[END_REF] showed how to improve on the order of these approximations for asymptotically normal linear statistics. For smooth non-linear statistics with non-zero asymptotic covariance, asymptotic normality still holds and Ith order approximations to their distributions and quantiles can be obtained by expanding their cumulants as power series in n -1 then applying Edgeworth-Cornish-Fisher technology. To do this two calculi are developed analogous but simpler to that of von Mises. In the case, where the mean is not changing, the corresponding tests can be viewed as one-and two-sided confidence intervals for the mean.

Set X 0 = 0 and let X = X 1 , X 2 , • • • be independent random variables in R p from some distribution F (x) say, with mean µ, and finite moments. (If p = 1 we denote the rth cumulant by κ r and set σ 2 = κ 2 .) These observations can be considered as a random process which may at some point go "out of control", changing their distribution. We define the average process of the observations as

M n (t) = n -1 S [nt]
(1.1)

ACCEPTED MANUSCRIPT

for S 0 = 0, S i = i j=1 X j and 0 ≤ t ≤ 1, where [x] is the integral part of x. A change in mean can be tracked by monitoring the average process via some functional of it, say T (M n ), referred to as a cusum statistic. We denote the mean of the average process (1.1) by

m n (t) = E M n (t) = µ[nt]/n → m(t) = µt (1.2)
as n → ∞ for 0 ≤ t ≤ 1. Given a functional T , θ = T (M n ) can be thought of as an estimate of θ = T (m). In this way we can if desired use θ to provide an estimate of µ. We shall assume that θ satisfies the standard cumulant expansion given by

κ r ( θ) = ∞ i=r-1 a ri n -i (1.3) for r = 1, 2, • • •.
For univariate data the most common prospective and retrospective two-sided cusum statistics and functionals are

A n = A(M n ) = n max k=1 |S k -kµ|/n, (1.4)
B n = B(M n ) = n max k=1 |S k -k Xn |/n (1.5) for A(g) = sup [0,1] |g(t) -tµ| and B(g) = sup [0,1] |g(t) -tg(1)|, where Xn = M n (1) = S n /n. If σ is a consistent estimate of the standard deviation σ, then as n → ∞ σ -1 n 1/2 {M n (t) -tµ} L → W (t), σ -1 n 1/2 {M n (t) -tM n (1)} L → W 0 (t)
for W (t) a Wiener process and W 0 (t) = W (t) -tW (1) a Brownian Bridge. So,

σ -1 n 1/2 A n L → sup [0,1] |W (t)|, (1.6)
σ -1 n 1/2 B n L → sup [0,1] |W 0 (t)|. (1.7)
See [START_REF] Billingsley | Convergence of Probability Measures[END_REF] and [START_REF] Anderson | Asymptotic theory of certain 'goodness of fit' criteria based on stochastic processes[END_REF].

These asymptotic results are easily extended to functionals like sup{|g(t)

-tµ 0 | -b(t)}, that is to test H 0 : µ = µ 0 versus H 1 : µ = µ 0 by rejecting H 0 if M n (t) -tµ 0 crosses a given boundary ±b(t).
An alternative is to use the L 1 norm. For example, the one-sided test of

H 0 : µ = µ 0 versus H 1 : µ > µ 0 one can use 1 0 {g(t)-tµ 0 -b(t)}dt, or equivalently 1 0 g(t)
dt. This is an example of the statistics and functionals we consider here. They include T (M n ) for

T (g) = 1 0 w(g(t), t)dt or 1 0 w(g(t), t)dg(t) (1.8)
for w(x, t) : R p × [0, 1] → R q , as well as smooth functions of such functionals. These functionals have the advantage of being asymptotically normal, unlike (1.6) and (1.7), and of having distribution, density and quantiles given by their Edgeworth-Cornish-Fisher expansions. In contrast, expansions for the distribution, density and quantiles of the cusum statistics (1.4)-(1.7), are not available.

ACCEPTED MANUSCRIPT

The aim of this note is to extend the work of [START_REF] Withers | Accurate tests and intervals based on linear cusum statistics[END_REF] for nonlinear functionals. In this case, the cumulant expansions in (1.3) and hence the Edgeworth-Cornish-Fisher expansions (Section 2, Withers and Nadarajah, 2009) are obtained in terms of functional derivatives of T at m using the Euler-McLaurin expansion. Section 2 defines the functional derivatives needed. There are two types depending on whether integration is with respect to dt or dg(t). These are analogous to those of von [START_REF] Mises | On the asymptotic distribution of differentiable statistical functions[END_REF] as extended in [START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF][START_REF] Withers | Nonparametric confidence intervals for functions of several distributions[END_REF]. We show that the chain rule applies to such derivatives. Examples include (1.8) and smooth functions of such functionals. Section 3 gives the cumulants of T (M n ) needed to obtain its Edgeworth-Cornish-Fisher expansions (Section 2, Withers and Nadarajah, 2009). For example, as n → ∞,

n var(T (M n )) → σ 2 T m (t) 2 dt
for m(t) of (1.2) and T g (t) the first derivative of T (g) in the calculus of Section 2.

Functional Derivatives

Since M n (t) is not a distribution function, the calculus of von [START_REF] Mises | On the asymptotic distribution of differentiable statistical functions[END_REF] does not apply to functionals of it. In this section we develop two differential calculi which can be applied to such functionals. The first calculus is for functionals like 1 0 w(g(t), t)dt when g(0) = 0 and for f ( wdg) without the condition g(0) = 0. The second calculus covers functionals like 1 0 w(g(t), t)dt without the condition g(0) = 0.

First Calculus

Consider functions g, h : I → R, where I ⊂ R p , and a functional T (g) in R. Suppose that for real ǫ,

T (g + ǫh) = T (g) + ∞ r=1 T r (g, h)ǫ r /r!, (2.1)
where

T r (g, h) = I • • • I T g (t 1 , . . . , t r )dh(t 1 ) . . . dh(t r )
is assumed to exist. We call T g (t 1 , . . . , t r ) the rth order (functional) derivative of T (g) with respect to the first calculus, or simply the first rth order derivative. It is made unique by the condition

T g (t 1 , . . . , t r ) is symmetric in t 1 , . . . , t r . ( 2

.2)

Restrictions on h may be required in order that (2.1) hold, as the examples show. That is, for a given T the expansion for T (g n ) -T (g) may place restrictions on eligible g ng.

Note 2.1 Condition (2.1) can be weakened in two ways. Firstly, it need only hold in a neighbourhood of zero, say |ǫ| < δ. This is satisfied by ǫ = n -1/2 if n > δ -2 . One can then take, for example, g(t) = m(t) and h(t) = n 1/2 (M n (t)m(t)), so that the left hand side of

ACCEPTED MANUSCRIPT

(2.1) is T (M n ). Secondly, the infinite Taylor series (2.1) can be replaced by a partial series with a remainder term, analogous to the usual ordinary Taylor series for a real function of a real variable, given in Section 3.6.1 of [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF].

Like von Mises' calculus, these calculi cannot be applied to sup functionals like that of (1.4) and (1.5).

Set ∆ t (s) = I(t ≤ s) = p i=1 I(t i ≤ s i ).
(2.3)

Taking λ 1 , • • • , λ r in R and h(s) = r i=1 λ i ∆ t i (s), where ∆ t i (•) is defined by (2.3), we see that T g (t 1 , . . . , t r ) is the coefficient of λ 1 • • • λ r in T r (g, h)/r!, so that T g (t 1 , . . . , t r ) is the coefficient of λ 1 • • • λ r ǫ r in T (g + ǫh)/r!. So, the first derivative is S(g) = T g (t) = lim ǫ→0 {T (g + ǫ∆ t ) -T (g)} /ǫ, (2.4)
where ∆ t is defined by (2.3). Similarly, by the method of proof of Theorem 2.1 of [START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF], it follows that

T g (t 1 , . . . , t r+1 ) = S g (t r+1 ) for S(g) = T g (t 1 , . . . , t r ). (2.5)
That is, the derivative of the rth derivative is the r + 1st derivative, just as for ordinary partial derivatives. The situation is analogous to, but simpler to the theory of derivatives of functionals of distribution functions developed by Von [START_REF] Mises | On the asymptotic distribution of differentiable statistical functions[END_REF] and refined by [START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF]. Note that in general 1 0 T g (t 1 , . . . , t r )dt 1 = 0, whereas for the von Mises derivative T F (x 1 , . . . , x r )dF (x 1 ) = 0: see (2.3) of [START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF].

Example 2.1 Suppose that p = 1 and T (g) = 1 0 w(t)dg(t), where w(t) : [0, 1] → R and g : [0, 1] → R are chosen such that T (g) exists. Set w r = ( w r ) 1/r for 0 < r < ∞. Then T g (t) = w(t) and higher derivatives are zero. So, T m dm = µw 1 = 0 in general.

Example 2.2 Suppose that p = 1 and T (g) = 1 0 W (t)g(t)dt. Then one needs to impose the condition h(0) = 0 giving T g (t) = 1 t W (s)ds or the condition h(1) = 0 giving T g (t) = t 0 W (s)ds. The higher derivatives are zero. The condition h(0) = 0 restricts application to T (g n ) -T (g) with g n (0)g(0) = 0, a condition satisfied by (g n , g) = (M n , m) of (1.1), (1.2). The condition h(1) = 0 restricts application to T (g n ) -T (g) with g n (1)g(1) = 0, a condition not satisfied by (g n , g) = (M n , m).

For the next examples it is helpful to use ∧ and ∨ to denote min and max:

t r∧ = r min i=1 t i , t r∨ = r max i=1 t i . ( 2
T g (t 1 , . . . , t r ) = w .r (g(s))∆ t 1 (s) • • • ∆ tr (s) = w .r (g(s))∆ tr∨ (s),

ACCEPTED MANUSCRIPT

where ∆ t 1 (•) is defined by (2.3). If h(s+) = h(s) and T (g) = w(g(s)g(1)) then

T g (t 1 , . . . , t r ) = w .r (g(s) -g(1))(∆ t 1 (s) -1) • • • (∆ tr (s) -1)
= (-1) r w .r (g(s)g(1))I(s < t r∧ ).

Example 2.4 Let ν be a measure on some space L. For s in L let T s (g) be a functional with finite rth derivative T sg (t If g is right-continuous (for example, g = M nm) and

S(g) = 1 0 w(g(s) -g(1), s)ds (2.8)
then S g (t 1 , . . . , t r ) = (-1) r tr∧ 0 w .r0 (g(s)g(1), s)ds.

(2.9)

For example, alternatives to the statistics A n and B n of (1.4) and (1.5) are

C n = C n (M n ) = n -3 n j=1 (S j -jµ) 2 w(j/n)/2, D n = D n (M n ) = n -3 n j=1 (S j -jS n /n) 2 w(j/n)/2 for C n (g) = 1 0 (g(s) -µ[ns]/n) 2 w([ns]/n)ds/2, D n (g) = 1 0 (g(s) -g(1)[ns]/n) 2 w([ns]/n)ds/2.
Their non-zero derivatives are

C ng (t 1 ) = 1 t 1 (g(s) -µ[ns]/n)w([ns]/n)ds, C ng (t 1 , t 2 ) = 1 t m2 w([ns]/n)ds, D ng (t 1 ) = 1 0 (g(s) -g(1)[ns]/n)w([ns]/n)e n (s, t 1 )ds, D ng (t 1 , t 2 ) =
C ′ n = C(M n ), D ′ n = D(M n ). If w(s) = 1 then one may show that C n -C ′ n and D n -D ′ n are O p (n -3/2 ) in probability.
Since D(g) does not have the form (2.8), (2.9) does not apply. This can be fixed by applying the method not to M n (t) but to

N n (t) = [nt] -1 [nt] i=1 X i ,
interpreted as 0 for t ≤ 1/n. However, κ( w 1 , • • • , w r ) for w i = N n (t i ) is more complicated than for w i = M n (t i ), so we do not pursue this idea beyond the next example.

Example 2.6 Consider

D ′ n = D ′ n (N n ) = n -1 n-1 j=1 (S j /j -S n /n) 2 w(j/n)/2,
where

D ′ n (g) = 1 n -1 (g(s) -g(1)) 2 w([ns]/n)ds/2.
The non-zero derivatives are

D ′ ng (t 1 ) = - 1 n -1 (g(s) -g(1))w([ns]/n)I(s ≤ t)ds, D ′ ng (t 1 , t 2 ) = 1 n -1 w([ns]/n)I(s < t m2 )ds. Note that D n [w] = D n and D ′ n [w] = D ′ n are equivalent in that D ′ n [w ′ ] = D n [w] for w ′ (t) = t 2 w(t). Note that 2n 3 D n [1] = 2n 3 D ′ n [t 2 ] = n j=1 (S j -j Xn ) 2
provides a simple alternative to

nB n = n max j=1 |S j -j Xn |,
since they provide L 2 and L ∞ norms of S jj Xn .

ACCEPTED MANUSCRIPT

Chain Rule

Since T g (t 1 , . . . , t r ) is an exact analogue of the ordinary rth order partial derivative, it satisfies the chain rule. We now spell this out in detail. By Comtet (1974, page 137), Faa di Bruno's chain rule for the derivatives of a function of a function is most simply stated in terms of the exponential Bell polynomials tabled on page 307 of [START_REF] Comtet | Advanced Combinatorics[END_REF]: for r ≥ 1, the rth derivative of a real function of a real function, say S(g) = f (T (g)), for g in R, T : R → R, f : R → R is given by

S (r) (g) = r k=1 f k B rk (T) (2.10) at f k = f (k) (T (g)) and T = (T 1 , T 2 , • • •)
, where T r = T (r) (g), the rth derivative of T (g). This extends to g and T (g) multi-dimensional by a simple reinterpretation. For example, for r = 4, (2.10) gives

S (4) (g) = f 1 T 4 + f 2 (4T 1 T 3 + 3T 2 2 ) + f 3 (6T 2 1 T 2 ) + f 4 T 4 1 .
For g in p dimensions and T (g) in q dimensions this becomes

S .i 1 i 2 i 3 i 4 (g) = f .j 1 T j 1 .i 1 i 2 i 3 i 4 + f .j 1 j 2 ( 4 i T j 1 .i 1 T j 2 .i 2 i 3 i 4 + 3 i T j 1 .i 1 i 2 T j 2 .i 3 i 4 ) +f .j 1 j 2 j 3 6 i T j 1 .i 1 T j 2 .i 2 T j 3 .i 3 i 4 + f .j 1 j 2 j 3 j 4 T j 1 .i 1 T j 2 .i 2 T j 3 .i 3 T j 4 .i 4 , (2.11)
where, for i

1 , i 2 , • • • in 1, 2, • • • , p S .i 1 i 2 ••• (g) = ∂ i 1 ∂ i 2 • • • S(g) for ∂ i = ∂/∂ g i .
Similarly, for the partial derivatives of T j (g) for j = 1, • • • , q; and for

j 1 , j 2 , • • • in 1, 2, • • • , q f .j 1 j 2 ••• (T (g)) = ∂ j 1 ∂ j 2 • • • f (y) at y = T (g) for ∂ j = ∂/∂ y j .
The convention is that j 1 , j 2 , • • •, the repeated pairs of suffixes in (2.11), are implicitly summed over their range 1, 2, • • • , q. The expressions N i in (2.11) mean summation over all N permutations of the i ′ s giving different terms. For example,

4 i T j 1 .i 1 T j 2 .i 2 i 3 i 4 = T j 1 .i 1 T j 2 .i 2 i 3 i 4 + T j 1 .i 2 T j 2 .i 1 i 3 i 4 + T j 1 .i 3 T j 2 .i 1 i 2 i 4 + T j 1 .i 4 T j 2 .i 1 i 2 i 3 .
This chain rule in turn has a simple reinterpretation for functions of functionals. Consider the real functional S(g) = f (T (g)), where T (g) is a q-dimensional functional on functions g : I → R, where I ⊂ R p , and f (T ) is an ordinary function from R q to R. Then, for example, the fourth order functional derivatives of S(g) are given in terms of the functional derivatives of T (g) by (2.11)

with S .i 1 i 2 ••• (g) replaced by S g (t 1 t 2 • • •) and T j.i 1 i 2 ••• replaced by T jg (t 1 t 2 • • •), the functional derivatives of T j (g): S g (t 1 t 2 t 3 t 4 ) = f .j 1 T j 1 g (t 1 t 2 t 3 t 4 ) + f .j 1 j 2 4 t T j 1 g (t 1 )T j 2 g (t 2 t 3 t 4 ) + 3 t T j 1 g (t 1 t 2 )T j 2 g (t 3 t 4 )) + f .j 1 j 2 j 3 6 t T j 1 g (t 1 )T j 2 g (t 2 )T j 3 g (t 3 t 4 )
+f .j 1 j 2 j 3 j 4 T j 1 g (t 1 )T j 2 g (t 2 )T j 3 g (t 3 )T j 4 g (t 4 ).

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Extension to Vector g

We now extend the results of this section to g, h, ǫ s-vectors. Let ǫh denote the s-vector

with ith component ǫ i h i for i = 1, • • • , s. For r in N s , where N = {0, 1, 2, • • •}, set |r| = s i=1 r i , ǫ r /r! = s i=1 ǫ r i /r i !.
Now consider functionals T (g) such that for all s-vector functions g, h on I ⊂ R p ,

T (g + ǫh) = r∈N s T r (g, h)ǫ r /r! for T r (g, h) = |r| I T g ({u ij , i = 1, • • • , s; j = 1, • • • r i }) s i=1 r i j=1 dh i (u ij ),
where |r| I denotes |r| integrals over I. The rth derivative, the above integrand, is again made unique by the condition that it is symmetric in each of its s rows. Taking

h i (u) = r i j=1 λ ij ∆ t i (u), where ∆ t i (•) is defined by (2.3), we see that this derivative is the coefficient of s i=1 r i j=1 λ ij in T r (g, h)/r!.
Let e is be the ith unit vector in R s . Then as in (2.4), for t in R the s first derivatives are

S i (g) = lim ǫ→0 {T (g + ǫe is ∆ t ) -T (g)} /ǫ,
where ∆ t is defined by (2.3). Similarly, (2.5) holds if t 1 , . . . , t r+1 and t 1 , . . . , t r are interpreted as ith rows with the other rows implicit and the same on both sides, that is Example 2.7 Take I ⊂ R p , g : I → R, w : I → R q , T (g) = wdg, S(g) = f ( wdg).

Then S g (t 1 , • • • , t r ) = f j1 •••jr ( wdg)w j 1 (t 1 ) • • • w jr (t r ).

Second Calculus

Consider functions g, h : I → R, where I ⊂ R p , and a functional T (g) in R. Suppose that for real ǫ, (2.1) holds, where now

T r (g, h) = I • • • I T g (t 1 , . . . , t r )h(t 1 )dt 1 . . . h(t r )dt r
is assumed to exist. We call T g (t 1 , . . . , t r ) the rth order (functional) derivative of T (g) with respect to the second calculus, or simply the second rth order derivative. It is made unique by the condition (2.2).

Let δ(t) denote the Dirac function on R. For t, s ∈ R q , set δ(t) = q i=1 δ(t i ) and δ t (s) = δ(ts). Taking λ 1 , • • • , λ r in R and h(s) = r i=1 λ i δ t i (s), we see that T g (t 1 , . . . , t r ) is the coefficient of λ 1 • • • λ r in T r (g, h)/r!. So, the first derivative is (2.4) with ∆ t replaced by δ t . Similarly, (2.5) and the chain rule hold.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Example 2.8 Take p = 1. Suppose that T (g) = 1 0 v(t)g(t)dt = vg say. Then T g (t) = v(t) and higher derivatives are zero.

Example 2.9 Take p = 1. Let s be any point in [0, 1]. Consider functions g : [0, 1] → R, and w : R → R. If T (g) = w(g(s)) then

T g (t 1 , . . . , t r ) = w .r (g(s))δ t 1 (s) • • • δ tr (s).
If T (g) = w(g(s)g( 1)) then

T g (t 1 , . . . , t r ) = w .r (g(s) -g(1))(δ t 1 (s) -δ t 1 (1)) • • • (δ tr (s) -δ tr (1)).
I ⊂ R p , g : I → R, v : I → R q , T (g) = I vg, S(g) = f ( I vg). Then S g (t 1 , • • • , t r ) = f j1 •••jr ( I vg)v j 1 (t 1 ) • • • v jr (t r ).
We note in passing that our method avoids the complexities of Gateaux, Fréchet and Hadamard derivatives. For more on these see, for example, Dudley (1992), [START_REF] Bednarski | Statistical expansions and locally uniform Fréchet differentiability[END_REF], Ren andSen (1991, 1995) and Pons and de Turckheim (1991).

Cumulant Expansions

Here we apply the first calculus of the previous section with p = 1, I = [0, 1] in order to obtain the leading coefficients {a ri } in the cumulant expansion (1.3) for θ = T (M n ). The higher order tests and confidence intervals can then be obtained as in Withers and Nadarajah (Section 2, 2009). By (2.1)

T (M n ) = T (m) + ∞ r=1 1 0 • • • 1 0 T m (t 1 , . . . , t r ) d(M n (t 1 ) -m(t 1 )) . . . d(M n (t r ) -m(t r ))/r!.
Cumulant expansions for T (M n ) could be obtained from this, following the method in [START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF], where cumulant expansions for T (F n ) were obtained for F n the empirical distribution of a random sample from F . However, in [START_REF] Withers | Nonparametric confidence intervals for functions of several distributions[END_REF] we showed that these could be obtained much more easily by applying the cumulant expansions of [START_REF] Withers | Second order inference for asymptotically normal random variables[END_REF] for t( w) by identifying w i with F n (x i ), i = 1, . . . , p and letting p → ∞. Here we take the same short cut with w i = M n (t i ). We begin with some notation.

For f : [0, 1] r → R, set (f ) r = 1 0 • • • 1 0 f (t 1 , • • • , t r )dt 1 • • • dt r . A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT for g 3 (t) = 1 0 T m (s)T m (s, t)ds.
In this way, discarding the terms in a rr of magnitude O(n -1 ) or less, one obtains

a 10 = T (m), a 21 = κ 2 [1 2 ], (3.1)
a 11 = κ 1 α 1 (T m ) + κ 2 [11]/2, (3.2)
a 32 = κ 3 [1 3 ] + 3κ {min(s 0 , s i )f 1 (s i )ds i }.

. 6 )

 6 If p > 1 then min and max are interpreted element-wise. Example 2.3 Take p = 1. Let s be any point in [0, 1]. Consider functions g : [0, 1] → R, and a smooth function w : R → R. Denote its rth derivative by w .r . If h(0) = 0 and T (g) = w(g(s)) then

1 0e

 1 n (s, t 1 )e n (s, t 2 )ds,A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPTwhere e n (s, t) = I(t ≤ s) -[ns]/n. Now replace C n (g), D n (g) by their limits s)g(1)s) 2 w(s)ds/2, and set

  Since the analog of Example 3.4 in Withers and Nadarajah (2009) holds for the second calculus, we have in the notation of Example 2.5, Example 2.10 Take

4 2 [ 1 , 2 , 3 , 1 0 1 0 1 0Example 3 . 1

 212311131 [START_REF] Withers | Accurate tests and intervals based on linear cusum statistics[END_REF] 234],where d 3 of a 33 is the bias term, obtainable from 3 ∆ abc of Withers (1983, page 67) asd 3 = α 1 (g) = (g(1)g(0))/2 for g(s) = 3κ 1 κ 3 T m (s, t)T m (t) 2 dt + 3κ 1 κ 2 2 {T m (s, t, u)T m (t)T m (u) +2T m (s, t)T m (t, u)T m (u)}dtdu. Suppose that T (g) = 10 w(g(s), s)ds. By (2.7) the rth derivative at g = m is 1 tr∨ f r (s)ds, where f r (s) = w .r0 (m(s), s). So, a 21 , a 11 , a 32 are given by (3.1)-(3.3) in terms of[1 r ] = • • • s r∧ r i=1 {f 1 (s i )ds i }, [11] = sf 2 (s)ds, α 1 (T m ) =f 1 (s)ds/2,[1, 12, 2] = f 2 (s 0 )ds 0 2 i=1

  1 , • • • , t r ). For S(g) = T s (g)dν(s),S g (t 1 , . . . , t r ) = T sg (t 1 , . . . , t r )dν(s).

	Putting these two examples together, we have
	1 tr∨	w .r0 (g(s), s)ds.	(2.7)

Example 2.5 Take p = 1. Let w(x, t) : R × [0, 1] → R have finite partial derivatives w .r0 (x, t) = (∂/∂ x ) r w(x, t). If g(0) = 0 and S(g) = 1 0 w(g(s), s)ds then S g (t 1 , . . . , t r ) =

  2 2[START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF][START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF][START_REF] Anderson | Asymptotic theory of certain 'goodness of fit' criteria based on stochastic processes[END_REF],(3.3)a 22 = κ 2 α 1 (T 2 m ) + κ 3 [1, 11] + κ 2 2 [12 2 ]/2 + κ 2 2 [1, 122] +2κ 1 κ 2 α 1 (g 3 ), a 43 = κ 4 [1 4 ] + 12κ 2 κ 3 [1 2 , 12, 2] + 12κ 3 2 [1, 12, 23, 3] +4κ 3 2 [1, 2, 3, 123].Note the formal equivalence, apart from the bias terms α k , with the expressions for the cumulant coefficients of T (F n ) given by Withers (1983, page 580). (The negative term in a 43 does not appear here due to the different form there for the rth order cumulants of w.) So, we obtain from[START_REF] Withers | Expansions for the distribution and quantiles of a regular functional of the empirical distribution with applications to nonparametric confidence intervals[END_REF] the coefficients needed to extend inference to O(n -2 ):a 12 = κ 1 α 2 (T m ) + κ 3 [111]/6 + κ 2 2 [1122]/8, a 33 = 3κ 4 [1 2 , 11]/2 + 3κ 2 κ 3 [1, 12, 22] + 3κ 2 κ 3 [1, 12 2 ] + 3κ 2 κ 3 [1, 2, 122] +3κ 2 κ 3 [1 2 , 122]/2 + 3κ3 2 [1, 2, 1233]/2 + κ 3 2 [12, 23, 31] +3κ 3 2 [1, 12, 233] + 3κ 3 2 [1, 23, 123] + d 3 , a 54 = κ 5 [1 5 ] + 20κ 2 κ 4 [1, 2 3 , 12] + 60κ 2 2 κ 3 [1, 2, 3, 12, 23] +60κ 2 2 κ 3 [1, 2 2 , 13, 23] + 30κ 2 2 κ 3 [1, 2, 3 2 , 123]

	+5κ 4 2 [1, 2, 3, 4, 1234] + 60κ 4 2 [1, 2, 34, 13, 24]
	+60κ

Acknowledgments

The authors would like to thank the Editor-in-Chief and the referee for carefully reading the paper and for their comments which greatly improved the paper.

ACCEPTED MANUSCRIPT

For π 1 , π 2 , • • • sets with elements in {1, 2,

3) holds with a ri given by [START_REF] Withers | Second order inference for asymptotically normal random variables[END_REF]

replaced by and t ij... replaced by T m (t 1 t 2 • • •). For example,

and the second term in a 11 is ij

The first term in a 11 is

for α i (g) of the Euler-McLaurin expansion (A.1) of Withers and Nadarajah (2009, Appendix A), since this implies that

where

By the corrigendum to [START_REF] Withers | Second order inference for asymptotically normal random variables[END_REF], the fifth term in a 22 is

= κ 1 κ 2 (g 3 (1)g 3 (0)) + O(n -1 )

ACCEPTED MANUSCRIPT

In particular, a 21 = κ 2 ds 1 ds 2 f 1 (s 1 )f 1 (s 2 )s 2∧ , 2a 11 = -µ f 1 (s)ds + κ 2 sf 2 (s)ds.