O Ertz

E Bocher
email: erwan.bocher@ec-nantes.fr

M Laurent

G Petit

D Rappo

STANDARD-CENTRIC AUTHORING AND PUBLICATION FOR CARTOGRAPHIC CONTENT

Keywords: standards, symbology encoding, map context, cartographic document, polypublishing/crossmedia

Within the nascent GeoWeb, democratization of GIS and cartographic practices is on the way with more and more map makers and map consumers looking for answers to their geographic and cartographic needs. This paper explains how standards are used as common languages to open the doors of collaborative edition of cartographic content. As more and more users get stuck in the making of maps, this paper describes some solutions to create adequate mapping tools adapted to various user profiles of skills. Especially for thematic mapping, until now the field of specialists, map makers want to take benefits of the analytical power of such maps. This paper brings some elements to anticipate this phenomena. Also, this paper explains how structuring cartographic content helps to make publication easier. The aim of this paper is to develop an innovative method of cartographic publication based on a polypublishing/crossmedia approach : one source, many media, many targets. It defines a cartographic document that needs to be adapted to these different media and targets while benefiting from all the winning cards of each.

INTRODUCTION

While the concept of Digital Earth [START_REF] Gore | The Digital Earth: Understanding our planet in the 21st Century[END_REF] is becoming reality through the nascent GeoWeb, democratization of GIS and cartographic practices is on the way with more and more map makers and map consumers looking for answers to their geographic and cartographic needs. Geospatial data are now widely available through web services and are the base material to create cartographic products useful for stakeholders. Based on best practices like use of standards and content structuring, the web is the ultimate platform to discover, explore, share, and create content by combining distributed informations and services. GeoWeb becomes a new paradigm for GIS and cartography to produce collective geospatial knowledge. By way of proof the growing of webmapping applications through geographic mashups (Programmableweb).

These applications offer new ways to explore and use geodata through the web medium, but printed maps stays important, while maps on mobile combined with location bases services becomes usual, and with no doubt, new channels for geodata and maps will appear in the future [START_REF] Crawford | Flexible Flat Panel Displays[END_REF]. In this context, it is a matter of fact that web mapping applications are tuned for one media and usually one target, for instance well done for screen display but poorly built up to provide also high quality printing. Still too often, for such applications, printing is secondary and when available, does not really fit the result expected, in term of resolution quality and layout. Many cartographic publication platforms are following the right direction with real printing functionalities (MapFish Print; GeoAdmin), but printing maps is still rather a desktop application functionality. However, online and offline worlds are coming closer so as cartographic publication platforms must now be able to adapt their contents in a multi-channels and multi-targets context. This is the main focus of this paper.

The red dot fever continues on one's way, leading to crowdsourcing, geospatial data infrastructures are built up, from local to global, and the next step is now to feed the geospatial knowledge infrastructure. It's time to create and share geospatial cognition, in the form of "mapping products" with the intend to synthesize informations, to communicate a message or to uphold a line of argument for taking decisions.

Thematic mapping represents an essential mechanism to absorb geodata whose volume is growing day by day on the GeoWeb. Even if this type of maps are more complex to build, until now the field of specialists, more and more map makers wants to take benefits of the analytical power of such maps. This paper brings some elements to anticipate this phenomena.

GEOSPATIAL KNOWLEDGE PRODUCTION

By extension of Spatial Data Infrastructure (Nebert, D., 2004), which concerns the discovering and sharing of geodata, a map product is a form of geospatial knowledge build up from geodata. And before mapping, a process can also produce the indicators and results to map. Finally, maps are published within cartographic documents.

The process to transform geodata into cartographic document for publication involves different types of actors :

1. data owner : to acquire and organize data according to a business logic;

2. data analyst : to process and transform "raw data" so as to extract indicators and results ready for mapping;

3. cartographer : to build map in respect of cartographic rules and best practices [START_REF] Bucher | The need for Web legend services[END_REF]; 4. cartopublisher : to transform maps into a cartographic document combined with other maps and enriched with external contents.

A geospatial knowledge infrastructure is driven by a collaborative context in which each actor pushes results of know-how in order to share with other actors :

1. push geodata on a GeoServer and push geo-statistical data on a StatServer;

2. push complex processes on a ProcessingServer;

3. push maps for rendering on a MapServer; 4. push cartographic document for publication on a CartoServer.

Within such an infrastructure, use of standards defined by Open Geospatial Consortium (OGC) is essential [START_REF] Iosifescu-Enescu | Cartographic Web Services And Cartographic Rules -A New Approach For Web Cartography[END_REF]. First, as shown in Figure 1, actors interact together through server components, so as standards are the engine oils for smooth-running transformation of geodata into cartographic document. Second, because actors want to work with their prefered tools, these have to be conform with these standards.

Figure 1. GeoCognition, many skills are in action to transform geodata into cartographic document While geodata can be accessed and processed in a standardized way using OGC standards like Web Feature Service (WFS), Geographic Markup Language (GML) and Web Processing Service (WPS), the following of this chapter will likely focus on mapping and publication of cartographic document using standards.

COLLABORATIVE CARTOGRAPHIC EDITION

In such a collaborative context, the use of a common language for cartographic description opens the doors of collaborative cartographic edition.

Map layer description

The OGC standard Symbology Encoding (SE) gives to several cartographers the ability to share the cartographic description of a map layer (see Figure 2).

Figure 2. Symbology Encoding, a common language to describe a map layer Indeed, the task of mapping indicators can involve both data analyst and cartographer profiles of actors. As an example, Figure 2 shows the use of SE, first to share symbologies that an expert have pre-defined, and to share a map layer built by an author but moderated (controlled) by an expert. As soon as the description is edited and validated, it is pushed for the rendering engine (FPS stands for Feature Portrayal Service) so as to style features retrieved from Web Feature Server.

Figure 3 shows a short diagram of Symbology Encoding model : a FeatureTypeStyle defines a pass to be applied on a map layer. Several Rules can then filter the features of the layer according spatial, logical, or comparison operators (e.g. st_area(the_geom) < 1000 AND name = 'Basel-Stadt' ). Each Rule is activated according domain constraints (e.g. 1:5000 < scale < 1:25000 ) and describes the symbolizer used to style the features that run in the filter. With OGC conform tools, a cartographer can then export the result of his(her) work according the Symbology Encoding standard model which has the capacity to formalize several types of maps. Currently, OGC is preparing the next release of the standard by reviewing Change Requests and Engineering Reports (OGC Charter document, 2010) that proposes improvements and new abilities. One Change Request is especially in interest for this chapter and concerns improvements for thematic mapping. To summarize, it suggests to add classification methods and a new type of graphic providing an internal way to describe charts like pie chart. In opposition with external approaches which does not guarantee a common language of description [START_REF] Mays | Using Styled Layer Descriptor definitions to display charts in a deegree Web Map Service[END_REF], the latest will extend the standard so as to be self-contained for such type of thematic maps.

Map context description

While the cartographer is defining the symbology of each layer, it is done with a map context in mind. Thus, a Map Context is a list of map layers that results in a map when put together. It is the base input material for the rendering phase to deliver a mapview. Rendering ability is either server-side or client-side.

Fro the image case, a layer can be a ready-to-view image (see Figure 4) produced by a server-side rendering engine with a Web Map Service (WMS) : a WMS request using default server-defined, a WMS-SLD * request with a custom symbology description by use of SE, or a FPS request combining a custom symbology with a WFS request (see Figure 2). The result is then an image in form of a formatted bitmap graphic like PNG or a vector graphic like SVG.

For the feature case, at the other side (Figure 4), a layer can be a WFS or TJS ( In a middle situation, a WPS request can provide a layer in the form of SLD description which combines inline features with a ready-to-use symbology description. Either it can then be rendered server-side through a FPS request or client-side. As an example, wpsPAL [START_REF] Ertz | Towards an automated label placement service[END_REF]) is such an OGC standards compliant cartographic service, it gives access to advanced labelling functionalities for any application that needs to produce cartographic content of high quality.

The second main important purpose of using a map context structured document concerns the description of metainformation useful for the final publications. It is essential to have title, source, date and so on for each map layer, but also layer-independent information such as the initial map extent, spatial reference system and some information describing the document itself such as title, abstract, source, author, etc.

The OGC Map Context standard is currently a candidate able to fulfil the description needs listed above. Thus several cartographers can share their work using this common language to describe a map. As soon as it is edited and validated, the map is pushed forward to be inserted in a cartographic document. At this time, it's up to the cartopublisher to complete the publication process.

* SLD stands for the OGC standard Styled Layer Descriptor

Authoring tools for map makers.

While more and more users get stuck in the making of maps, it's essential to create adequate mapping tools adapted to various user profiles of skills. Especially for thematic mapping which implies the respect of precise cartographic rules, some wizard tools would guide the creation of such type of maps but could also help to control the aptness of the map expression.

For thematic maps, with the OGC Geographic Data Attribute Set formatting standard (GDAS), used in conjunction with TJS, editing tools are able to get the attribute data but also its associated meta-data. This attribute typology (nominal, ordinal, count, measure) facilitates then automated determination of appropriate forms of mapping.

Also, some smart interactive visual widgets could help map makers. For instance concerning the creation of a choropleth map, the systematic use of an histogram to visualize the distribution is first required. Moreover, when choosing a classification method, providing a way to adjust threshold and visualize "live on" the results on the map would enrich the user experience and ease the work. Also, a systematic use of a "ColorBrewer"-like tool [START_REF] Brewer | ColorBrewer. The Pennsylvania State University[END_REF]) helps for selecting good colour schemes. And another smart mapping widget could help to automate and assist the user for the calculation of proportional symbol sizes.

CARTOGRAPHIC POLYPUBLISHING/CROSSMEDIA

Using common languages for map layer symbology and map context descriptions allows to put in practice collaborative edition of cartographic content. But also, structuring content make publication easier. This section aims to develop an innovative method of cartographic publication based on a polypublishing/crossmedia approach : one source, many media, many targets.

Cartographic document definition

A map context describes the information needed to create a map while a cartographic document can encapsulate several maps. In a multiple maps situation, there are then some documentspecific information like document title, sub-title, author, date, logo and so on. Moreover, a cartographic document is also the right place to enrich encapsulated maps with external content like text, images, etc. As maps are visualizations to support a line of argument expressed by the cartographic document, often some text is joined to bring some explanations. Below are use cases to illustrate this definition.

CartoExplorer case : considering the data owner has pushed a GPS track with its waypoints and some descriptive attributes (name, author, distance, height variation, etc) on a GeoServer.

The cartographer build up the map context :

•

define one or more possible base layers;

• define one overlay for the track and one for waypoints;

• define a layer of point of interests (POI) provided by a remote POI search service through a WPS request;

• define the symbology of each layer by use of Symbology Encoding according eventually a set of pre-defined symbologies;

• define initial map extent and spatial reference system;

• describe meta-information, layer name, date, author, source, etc.

The cartopublisher build up the cartographic document :

• insert one or more prepared map contexts;

• enrich the document with photos in relation with the cartographic theme;

• describe more meta-information if needed, especially if multiple maps are encapsulated.

Thematic mapping case : considering data owners have pushed geodata on a GeoServer and some data statistics on a StatServer with common geocodes.

The cartographer build up the map context :

• define a base layer;

• define thematic overlays through TJS requests;

• define the thematic symbology by use of Symbology Encoding with thematic mapping level of conformance;

• define initial map extent and spatial reference system.

The cartopublisher build up the cartographic document :

• insert one or more prepared map contexts;

• enrich the document with analytical explanations, article from newspaper, graphical charts, etc.;

• describe document-specific meta-information.

In collaboration with the cartographer, the cartopublisher creates the cartographic logic. As soon as the cartographic document is edited and validated, it represents one source ready for publication.

Polypublishing/crossmedia approach

This approach is synonym to multi-channels and multi-targets publications, meaning that a source document is described one time, containing all the base material, and published for many media and many targets. Crossmedia means for print medium, for screen display, for mobile, etc. each device-specific.

Polypublishing means from one source to many targets depending on final user needs and publication context.

Let's apply this approach for cartographic content publication.

For different channels :

• print : high-quality maps and statical, offline content;

• screen : online and interactive content through a webmapping application.

For different targets :

• only maps : rapid and synthetic for publication in a newspaper;

• maps with tables and charts : detailed and analytical for stakeholders.

Media and targets have their specificities, advantages and drawbacks. Print medium offers high-quality of rendering which enforces the principle to describe a cartographic document in complete independence with the final publication device. This means that all measures, especially for symbology description have to be defined with absolute units (e.g. not pixel). At the other side, web medium offers low resolution rendering but have a strong advantage with the ability to offer interactions for rich user experience when exploring maps.

Thus, a cartographic document needs then to be adapted to these different media and targets, and it is to benefit from all the winning cards of each.

Publishing process

A cartographic document must be transformed to be adapted according media and targets. The principle is that such a transformation knows the advantages and drawbacks of the targeted media so as to adapt the publication the best according specificities like resolution, interactive context, online/offline content, etc.

It's up to the cartopublisher to configure these transformations which first need a formatting language to fit the final targeted container. For print medium, the final container is a printer, thus a page description language is used like PDF, XPS or PostScript to formalize the result. For web medium, the final container is a screen display running a browser with a webmapping framework, using W3C standards and JavaScript language, like OpenLayers does.

Finally, each information element defined by the cartographic document has to take place within a presentation through a layout template. For instance, either it is a page layout to be printed or a web page layout to display on a screen. Using a template system allows then the reuse of the same template to publish other cartographic documents. It makes sense for an organization which has a pre-defined corporate template so as to standardize all cartographic publications produced by its members. Several technologies can be used to implement and execute the transformation chain. Table 1 shows two examples, one for print medium using the well-known XSL-FO language. It is a W3C recommendation that stands for Extensible Stylesheet Language Formatting Objects. It is a language for formatting XML data for output to screen, paper or other media. Mainly used for printed publications in conjunction with Apache FOP processor which produces output formats like PDF, PS, etc.

Templating

The second publishing process from Table 1 shows the transformation of a cartographic document towards a webmapping application in a client/server publication context. Thus, server-side, Smarty, a template system for PHP manages the presentation while HTML, CSS and Javascript are the used formatting languages for running within OpenLayers, a wellknown AJAX mapping framework that makes it easy to put a dynamic map in any web page.

Among other technologies, GeoPDF and KML would be also interesting candidate for publication.

Let's zoom into the process. For each map context described, it has to be drawn, and this is where it comes to know about the targeted resolution which is essential for high-quality printing.

The good news for a WMS map layer is that the OGC Standard Working Group is currently considering suggestions to have a resolution parameter (i.e. dpi) within a WMS request. This, in addition to absolute units for Symbology Encoding, would then lead to a device-independent context. Also, one other solution would be to request an image using a scalable format like SVG.

Next step is to apply the template and to place the different content elements at the right place using the right form. This step is driven by a publication layout combined with configurable publication modules. Together they define the final user interface.

Figure 5. different types of map-legend

Modules are containers able to present content in a given form depending on how they are configured. For instance, a maplegend module is in charge of generating map-legends according Symbology Encoding descriptions and configuration parameters like the type of representation (Figure 5).

It goes the same for the configurations of an orientation module to display the north arrow, a scale-bar module, etc. It's up to the cartopublisher to configure the publication through these modules.

Independently of the media or target, transformations are built so as to automate and guide the cartopublisher to apply the cartographic best practices, like the mandatory presentation of a title, a scale-bar, a north arrow, all map-legends, sources, dates, authors and so on.

One source, many specific transformations

As introduced previously, transformations adapt the content publication according final medium and target. While transformations will have a focus on rendering quality for print medium, one major advantage of a webmapping publication is the ability to provide map interactions. This means that modules can deploy specific elements of interactivity. The interactive context can be used to activate the usual navigation tools to zoom, to pan and to switch layers, but also some interactions to better analyse a thematic map. Figure 6 shows possible interactions between a map-view, map-legends, a data table and a graphic chart to help the user to identify all the information about an object selected on the map. The equivalent printed publication can present all this information, but in a statical manner.

To illustrate the method, let's come back to the first case described previously with its cartographic document and complete the publishing process with some transformations.

Example of transformations

The webmapping transformation needs a template layout to organize the presentation and distribute the content among the usual mapping modules.

First, let's consider a classical webmapping layout and some modules (Figure 7) :

1. Annotation module to style meta-information like title, date, author, etc.;

2. Zoom-bar module to allow map navigation;

3. Layer-switcher module to control layer visibility;

4. Feature-info module to allow object selection on the map and visualization of its descriptive attributes (needs configuration to associate a label to each attribute);

5. Map-legend module providing visual description of visible layers with sources information and interaction to zoom into a specific layer;

6. Photo strip module to provide an interactive way to explore external content like photos. One complementary module can be added for another target, an advanced publication, by including a profiler visualization (Figure 8). Interactive synchronization between the track on the map and the profile is a way to enrich user experience.

Figure 8. A map synchronized profiler module For a printed publication of the same content, there will be no interactivity, but the cartopublisher has to think the publication in a different way.

Figure 9. Multi-pages printed publication First, the map is printed on a full landscape page so as to intensify the reading experience. Then one page per layer is printed to lay down feature attributes. The latest will need a clean template, not just a raw table. For instance, if the layer is styled using a categorization, the table should list ordered features by classes with colours.

CONCLUSIONS

Use of standards for map layer symbology and map context descriptions is essential to put in practice collaborative edition of cartographic content. Moreover, cartographic editing tools can be improved to better suit user skills by a systematic use of attribute typology and smart mapping widgets to automate or at least guide the cartographic production.

Cartographic content structuring favour publication over the GeoWeb and allows easy syndication for derived geospatial analysis and knowledge production. Moreover, it opens the doors for multi-channels and multi-targets cartographic content publication. Such an approach is central for the design of a platform when it comes to publish cartographic content over various media, from usual high-quality printed maps to webmapping applications on mobile, and for various targets, depending on targeted technologies and final user profiles, from stakeholder to newspaper reader.

Figure 3 .

 3 Figure 3. Symbology Encoding model Several symbolizers can be defined : 1. AreaSymbolizer for the representation of areas given outline and filling values; 2. LineSymbolizer for the representation of lines given a type of stroke; 3. PointSymbolizer for the representation of points to put a graphic on (i.e. an image); 4. TextSymbolizer for the representation of text labels according placement options.

Figure

  Figure 4. Map Context

Figure 6 .

 6 Figure 6. Specific elements of interactivity for thematic mapping

Figure 7 .

 7 Figure 7. Webmapping layout and modules

Table 1 .

 1 Two publishing processes with technologies, one for print medium, one other for screen display

		Language	Engine	Result
	XSL	XSL-FO	FOP	PDF file
	Server-side	HTML/CSS	OpenLayers	Client-side
	templating	+JavaScript	AJAX mapping	webmapping
	(e.g. Smarty)		framework	application


Luc Pronzato
email: pronzato@i3s.unice.frrendas@i3s.unice.fr

Joaõ Rendas

Bayesian Local Kriging *

Keywords: prediction, interpolation, computer experiments, random field, non-stationary process

We consider the problem of constructing metamodels for computationally expensive simulation codes; that is, we construct interpolation/prediction of functions values (responses) from a finite collection of evaluations (observations). We use Gaussian process modeling and Kriging, and combine a Bayesian approach, based on a finite set of covariance functions, with the use of localized models, indexed by the point where the prediction is made. Our approach does not yield a single generative model for the unknown function, but by letting the weights of the different covariance functions depend on the prediction site, it gives enough flexibility for predictions to accommodate to non-stationarity. Contrary to Kriging prediction with plug-in parameter estimates, the resulting Bayesian predictor is constructed explicitly, without requiring any numerical optimization. It inherits the smoothness properties of the covariance functions that are used and its superiority over the plug-in Kriging predictor (sometimes also called empirical-best-linear-unbiased predictor) is illustrated on various examples, including the reconstruction of an oceanographic field from a small number of observations.

Introduction and motivation

The usual approach to Gaussian Process (GP) modeling and Kriging prediction raises two major issues: (i) stationarity is often a too strong assumption but seems hardly avoidable when a single realization of the random field is observed; (ii) the estimation of the kernel parameters that specify the correlation between distant observations is problematic, taking the uncertainty on these estimates into account in the construction of predictions is difficult and requires either heavy Monte-Carlo calculations or relies on (sometimes crude) approximations based on the asymptotic behavior of the estimators.

The classical plug-in approach, with which we shall compare, consists in predicting with a correlation function where unknown parameters are replaced by their estimated values, usually by Maximum Likelihood (ML). This method is simple but may produce poor results, in particular (a) when the modeled phenomenon is strongly non-stationary, (b) when an unlucky poor sampling wrongly suggests that the process is exaggeratedly smooth (this corresponds to the notion of deceptive function, see [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]), or (c) far away from the design points where the process is observed, because Kriging prediction is then given by the global trend and may be locally very inaccurate. Such situations will be considered in Section 3.

To address difficulties (i) and (ii) above, we propose a local-Kriging approach that combines two features.

First, the approach is Bayesian. It relies on a finite set of L GP models {Z (x)} L =1 , for x varying in a given set X , Z (•) having covariance function C . We then consider that f (•) is the sample path of a process Z s (•) such that s = with some probability w . Starting with prior weights w 0 (for instance uniform), we can update them into w n after n observations z n = (Z(x 1 ), . . . , Z(x n )) have been collected, and hence construct a Bayesian predictor x ∈ X -→ ηn (x) based on the L models. This prediction and its posterior squared error can be constructed explicitly when we assume a linear parametric trend g (x)β, with g(•) a known vector of functions (the usual framework for universal Kriging), and a hierarchical prior β|σ 2 ∼ N (β 0 , σ 2 V 0 ) and σ 2 ∼ inverse chi-squared, common to the L models. This is rather standard, see Santner et al. (2003, Chap. 4). Notice that, due to the dependence of the w n in z n , this Bayesian predictor depends non linearly in z n . The method differs from the approach based on mixture of kernels in [START_REF] Ginsbourger | Discrete mixture of kernels for krigingbased optimization[END_REF], which is not Bayesian, and it extends [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF] by allowing the possible assignment of prior weights to different covariance structures. Several examples (Section 3) will illustrate that a small number L (typically L ≈ 4) of isotropic covariances is enough to obtain satisfactory results, at least for the two-dimensional sets X considered.

The second feature of our construction is meant to account for non-stationarity. Instead of proposing a unique view of the process realization, or in other words, of the phenomenon to be modeled, we consider that observers at different locations t and t may contemplate different models. We thus condition all process characteristics, in particular the covariances C , by the location t where prediction is made. This yields L covariance functions C |t for each t, and a Bayesian predictor ηn (•|t). By construction, the prediction ηn (x 0 |t) at x 0 is only valid locally, for x 0 close to t, but the predictor t ∈ X -→ ηBLK,n (t) = ηn (t|t), which we call Bayesian Local Kriging (BLK), inherits continuity and interpolating properties from the properties of the covariances C |t . The method differs from previously proposed local Kriging approaches, see for instance [START_REF] Lam | A novel meshless approach -Local kriging (LoKriging) method with two-dimensional structural analysis[END_REF]; [START_REF] Sun | Analysis and prediction of soil properties using local regression-kriging[END_REF][START_REF] Nguyen-Tuong | Model learning with local Gaussian process regression[END_REF]. It does not yield a single generative model for the process, but by assigning different posterior weights w n (t) to different locations t, it gives enough flexibility for predictions to accommodate to nonstationarity.

Construction

The objective is to reconstruct/predict/interpolate a function f (

•): x ∈ X -→ f (x) ∈ R from a collection of observations f (x 1 ), . . . , f (x n ) at sites x 1 , . . . , x n in X , with X a compact subset of R d , d ≥ 1.
Real-valued GP models will be used to predict f (x 0 ) at some unsampled x 0 by Kriging. Since there is no reason to consider f (•) as the sample path of a stationary GP, we shall use collections of local models, each of them defining a local prior representation for f (•). We shall write

Z(•) ∼ GP(µ(•), σ 2 C(•, •)) when Z(•) is a GP satisfying E{Z(x)} = µ(x) , ∀x ∈ X , and E{[Z(x) -µ(x)][Z(x ) -µ(x )]} = σ 2 C(x, x ) , ∀(x, x ) ∈ X 2 .

Local models

For the sake of simplicity, for the moment we consider processes with zero mean; see the Appendix for the case where a linear parametric trend is present. At any site t ∈ X , we shall use a set of L local GP models GP(0, σ 2 C |t (•, •)) for f (•). Non stationarity is introduced by letting C |t (x, x ) depend on t, which makes the model local. We shall use covariance functions of the form

C |t (x, x ) = k 0, (x -x ) k 1/2 1, (x -t)k 1/2 1, (x -t) , (1)
where the k i, (•), i = 0, 1, = 1, . . . , L, are stationary kernel functions, continuous at 0 and such that k i, (0) = 1. The covariance

C |t (x, x ) satisfies C |t (x, x) = k -1 1, (x -t) and C |t (t, t) = 1
. By taking k 1, (x -t) decreasing in x -t , we obtain a prior for f (•) which may become quite vague for x far enough from t (however, the correlation corresponding to C |t (x, x ) is k 0, (x -x ) for any t and is thus stationary).

The Bayesian local Kriging predictor

Consider a given site t ∈ X . With each model GP(0,

σ 2 C |t (•, •)), = 1, . . . , L, with C |t (•, •)
given by (1), we associate a prior probability w 0 . That is, we consider that, seen from t, f (•) is the sample path of a stochastic process Z(•) such that

Z(•)|s(t), σ 2 , t ∼ GP(0, σ 2 C s(t)|t (•, •)) , with Prob{s(t) = } = w 0 , = 1, . . . , L , (2)
with some prior distribution on σ 2 (common to all components), having density ϕ 0 (•) with respect to the Lebesgue measure on R + . We thus have a finite mixture of GP's at each t, and different mixtures at t and t = t.

Remark 1 (Linear combination of GP) Note the difference between this Bayesian construction and the more usual consideration of a linear combination of processes, that is, [START_REF] Nott | Estimation of nonstationary spatial covariance structure[END_REF]. In that case, assuming that E{Z (x)Z (x )} = 0 for all (x, x ) ∈ X 2 when = , we obtain that E{Z(x)Z(x )} = σ 2 L =1 w (x)w (x )C (x, x ), which makes Z(•) non stationary also when C (x, x ) = C (x -x ) for all (x, x ). This model may seem simpler to handle than the Bayesian hierarchical one considered above. However, predictions require the estimation of the weight functions w (•), which is hardly possible when a single realization of Z(•) is available and no prior parametric model is available for w (x).

Z(x) = L =1 w (x) Z (x), where Z (•) ∼ GP(0, σ 2 C (•, •)), see, e.g.,
After n evaluations of f (•) at sites

ξ n = (x 1 , . . . , x n ) in X , the likelihood L (z n | , σ 2 , t) of z n = (f (x 1 ), . . . , f (x n )) for the model GP(0, σ 2 C |t (•, •)) is given by L (z n | , σ 2 , t) = 1 σ n (2π) n/2 det 1/2 K n ( |t) exp - 1 2σ 2 z n K -1 n ( |t)z n , where K n ( |t) is the matrix with elements {K n ( |t)} i,j = C |t (x i , x j ), i, j = 1, . . . , n. Notice that K n ( |t) = D n ( |t)K n,0 ( )D n ( |t), where {K n,0 ( )} i,j = k 0, (x i -x j ) and D n ( |t) = diag{k -1/2 1, (x i -t) , i = 1, . . . , n} , see (1)
. Therefore, as we stressed before, for suitable kernels k 1, (•), the information carried by observation z i = f (x i ) decreases with the distance from t to x i . Using Bayes rule, the marginal likelihood

L (z n | , t) = R + L (z n | , σ 2 , t) ϕ 0 (σ 2
) dσ 2 is used to compute the posterior weights

w n (t) = w 0 L (z n | , t) L =1 w 0 L (z n | , t) , = 1, . . . , L ; (3)
that is, the posterior distribution of s(t).

Remark 2 The different covariance functions C |t (•, •), = 1, . . . , L, may correspond to different parameter values θ l in a given parameterized kernel C |t (•, •|θ), typically different length scales. The method can then be extended straightforwardly to the infinite mixture case, with a continuous limit that gives some prior density π 0 (•) to θ. The posterior is then

π n (θ) = π 0 (θ)L (z n |θ, t) π 0 (θ)L (z n |θ, t) dθ .
However, in general the integral on the denominator cannot be calculated analytically, and one must resort to MCMC methods or use the Laplace approximation. Note that, besides being more easily tractable, the finite mixture model considered here allows one to let the different C |t (•, •) represent different correlation characteristics (isotropy, smoothness, etc.).

Let η 0 denote the prediction of the unobserved value Z(x 0 ) for the local model at site t. Its posterior squared prediction error is

PSPE(η 0 ) = E{[Z(x 0 ) -η 0 ] 2 |z n , t} ,
where the expectation is with respect to Z(x 0 ) given z n and t, with Z(x 0 )|s(t), σ 2 , z n , t ∼ GP(0, σ 2 C s(t)|t (•, •)|z n ), Prob{s(t) = } = w 0 , l = 1, . . . , L (which gives again a finite mixture of GP), and

σ 2 ∼ ϕ 0 (•). Denote ηn (x 0 |t) = E{Z(x 0 )|z n , t} . (4)
Then, PSPE(η 0 ) = var{Z(x 0 )|z n , t} + [η n (x 0 |t) -η 0 ] 2 , which is minimum for η 0 = ηn (x 0 |t). The associated posterior squared prediction error equals the posterior variance

PSPE(η n (x 0 |t)) = var{Z(x 0 )|z n , t} = L =1 w n (t) E{[Z(x 0 ) -ηn (x 0 |t)] 2 |z n , , t} = L =1 w n (t) E{[Z(x 0 ) -ηn (x 0 | , t)] 2 |z n , , t} + L =1 w n (t) [η n (x 0 | , t) -ηn (x 0 |t)] 2 ,
where we have denoted ηn

(x 0 | , t) = E{Z(x 0 )|z n , , t} . Since ηn (x 0 |t) = L =1 w n (t)η n (x 0 | , t), we obtain PSPE(η n (x 0 |t)) = L =1 w n (t) var{Z(x 0 )|z n , , t} + var{d n (x 0 |t)} , (5)
where d n (x 0 |t) denotes the discrete distribution that allocates weight w n (t) to ηn (x 0 | , t).

Due to the use of localized covariance functions, see (1), the prediction ηn (x 0 |t) only makes sense for x 0 close to t. We call Bayesian Local Kriging (BLK) the predictor t ∈ X -→ ηBLK,n (t) = ηn (t|t) .

Note that it depends nonlinearly on z n even in situations where each prediction ηn (t| , t) is linear in z n , since the weights w n depend on z n through the marginal likelihood L (z n | , t), see (3). Also note that ηn (t|t) inherits the smoothness properties of k 0, (•) and k 1, (•). The posterior squared prediction errors, or posterior variances, PSPE(η BLK,n (t)) = var{Z(t)|z n , t}, t ∈ X , can be used to measure the precision of predictions made after the collection of observations z n . The preposterior variances, or mean-squared prediction errors,

MSPE(η BLK,n (t)) = E{var{Z(t)|z n , t}|t} , t ∈ X , (6)
can be used for experimental design: one may choose a set of locations ξ n = (x 1 , . . . , x n ) ∈ X n that ensures a precise prediction of the behavior of f (•) over X by minimizing

Φ M (ξ n ) = max t∈X MSPE(η BLK,n (t)) or Φ I (ξ n ) = X MSPE(η BLK,n (t))ζ(dt) , (7)
with ζ(•) some suitable measure on X . As shown in Section 2.3, the predictor (4) and its prediction error (5) can be expressed explicitly when the prior density ϕ 0 (•) on σ 2 is suitably chosen (conjugate prior).

An inverse chi-square prior for the process variance

Suppose that Z(•) satisfies (2) with σ 2 having the inverse chi-squared distribution (common to all components),

ϕ 0 (σ 2 ) = ϕ σ 2 0 ,ν 0 (σ 2 ) = (σ 2 0 ν 0 /2) ν 0 /2 Γ(ν 0 /2) exp[-ν 0 σ 2 0 /(2σ 2 )] σ 2(1+ν 0 /2) , ( 8
) so that E{1/σ 2 } = 1/σ 2 0 , var{1/σ 2 } = 2/(ν 0 σ 4 0 ), E{σ 2 } = σ 2 0 ν 0 /(ν 0 -2) (for ν 0 > 2) and var{σ 2 } = 2σ 4 0 ν 2 0 /[(ν 0 -2) 2 (ν 0 -4)] (for ν 0 > 4). Direct calculations give ∞ 0 L (z n | , t)ϕ 0 (σ 2 ) dσ 2 = 1 (2π) n/2 det 1/2 K n ( |t) (σ 2 0 ν 0 /2) ν 0 /2 Γ(ν 0 /2) Γ(ν n /2) (σ 2 n| ,t ν n /2) νn/2 ,
where ν n = ν 0 + n and

σ 2 n| ,t = ν 0 σ 2 0 + nσ 2 n ( |t) ν 0 + n (9) with σ2 n ( |t) = z n K -1 n ( |t)z n /n
, the Maximum-Likelihood (ML) estimator of σ 2 given z n , , t. Also, given z n , and t, σ 2 has the inverse chi-squared distribution ϕ n| ,t (•) = ϕ σ 2 n| ,t ,νn (•), see (8).

The posterior mean ηn (x 0 | , t) corresponds to the ordinary-Kriging predictor, also called Best Linear Unbiased Predictor (BLUP), c n (x 0 | , t)z n for the model GP(0,

σ 2 C |t (•, •)|z n ), with, c n (x 0 | , t) = K -1 n ( |t)k n (x 0 , |t) (10)
and

{k n (x 0 , |t)} i = C |t (x 0 , x i ), i = 1, . . . , n. Therefore, ηn (x 0 |t) = L =1 w n (t)c n (x 0 | , t)z n . (11)
The variance var{Z(x 0 )|z n , , σ 2 , t} is the ordinary-Kriging variance for GP(0,

σ 2 C |t (•, •)|z n ), var{Z(x 0 )|z n , , σ 2 , t} = σ 2 ρ 2 n (x 0 | , t), with ρ 2 n (x 0 | , t) = C (x 0 , x 0 |t) -k n (x 0 , |t)K -1 n ( |t)k n (x 0 , |t) (12)
(note that it depends on the design

ξ n = (x 1 , . . . , x n ) but not on z n ). Since E{Z(x 0 )|z n , , σ 2 , t} = ηn (x 0 | , t) does noes not depend on σ 2 , var{Z(x 0 )|z n , , t} = E{σ 2 |z n , , t} ρ 2 n (x 0 | , t) = σ 2 n| ,t ν n /(ν n -2) ρ 2 n (x 0 | , t) .
From ( 5) and ( 11), the posterior squared error of the prediction ηn (x 0 |t) at x 0 equals

PSPE(η n (x 0 |t)) = L =1 w n (t) σ 2 n| ,t ν n /(ν n -2) ρ 2 n (x 0 | , t) + z n Ω n (x 0 |t)z n , (13)
where Ω n (x 0 |t) is the variance-covariance matrix Var{D n (x 0 |t)} with D n (x 0 |t) the discrete distribution that allocates weight w n (t) to the vector c n (x 0 | , t), = 1, . . . , L.

The preposterior variance (mean-squared prediction error) of BLK at t is (for ν 0 > 2)

MSPE(η BLK,n (t)) = L =1 E{w n (t) σ 2 n| ,t ν n /(ν n -2)} ρ 2 n (t| , t) + E{z n Ω n (t|t)z n |t} = σ 2 0 ν 0 /(ν 0 -2) L =1 w 0 (t) ρ 2 n (t| , t) + E{z n Ω n (t|t)z n |t} . ( 14
) Note that MSPE(η BLK,n (t)) = E{σ 2 } ρ 2 n (t|t) when L = 1. The term E{z n Ω n (t|t)z n |t} = L =1 w 0 E{var{d n (x 0 |t)}| , t} , see (5)
, plays a role similar to that of the correcting term added to the Kriging variance in [START_REF] Harville | Mean squared error of estimation or prediction under a general linear model[END_REF][START_REF] Zimmerman | Mean squared prediction error in the spatial linear model with estimated covariance parameters[END_REF][START_REF] Abt | Estimating the prediction mean squared error in gaussian stochastic processes with exponential correlation structure[END_REF][START_REF] Zhu | Spatial sampling design under the infill asymptotic framework[END_REF].

The construction of an experimental design optimal in terms of Φ M (•) or Φ I (•), see ( 7), would require the evaluation of E{z n Ω n (t|t)z n |t} and deserves further investigations. A sequential approach facilitates the construction: suppose that n observations have already been collected, with ηBLK,n (t) the prediction at t, see ( 11), and PSPE(η BLK,n (t)) the associated posterior squared prediction error, see ( 13); a reasonable choice then places next observation at x n+1 where PSPE(η BLK,n (x)) is maximum.

In this section we have considered GP with zero mean. The presence of a linear parametric trend (universal Kriging) is considered in Section 5: expressions (11), ( 13) and ( 14) remain valid, but with different values for ν n , σ 2 n| ,t , c n (x 0 | , t) and ρ 2 n (x 0 | , t).

Examples

We shall consider isotropic kernels k i, (•) of the form k i, (x -x ) = K γ,θ ( x -x ), i = 0, 1, where K γ,θ (τ ) denotes a Matérn covariance function, with

K γ,θ (τ ) = (θτ + 1) exp(-θτ ) for γ = 3/2 (θ 2 τ 2 /3 + θτ + 1) exp(-θτ ) for γ = 5/2 (15)
as special cases, see, e.g., Stein (1999, pp. 31, 48). For fixed γ and θ, K γ,θ (τ ) is a decreasing function of τ ∈ R + , with a correlation K γ,θ (τ ) = 20% for τ 2.9943/θ when γ = 3/2 and for τ 3.9141/θ when γ = 5/2. The design space X is often renormalized to [0, 1] d , so that prior guesses on reasonable inverse correlation lengths θ can be set depending on the assumed smoothness of the function f (•). One dimensional processes with covariance K γ,θ (τ ) are m times mean-square (and almost surely) differentiable if and only if γ > m, see Cramér and Leadbetter (1967, p. 185), Stein (1999, p. 32), and K 3/2,θ (•) (respectively K 5/2,θ (•)) yields one-time (resp. two-times) isotropic differentiable processes on R d .

Example 1: average performance. We simulate various realizations of stochastic processes on X = [0, 1] 2 , stationary or not, using a linear combination of processes as mentioned in Remark 1,

Z(x) = [1 -w(x)]Z 1 (x) + w(x)Z 2 (x) , (16)
where Z 1 (•) and Z 2 (•) are stationary, with mean, variance and covariance functions βi , σ 2 i and (15). For the stationary case we shall use w(x) ≡ 1, so that Z(x) = Z 2 (x) for all x; realizations of a non-stationary process will be generated with w(x) = w n-s (x) given by the product of two sigmoid functions,

K γ i ,θ i (•), i = 1, 2, see
w n-s (x) = exp[a(x 1 -1/2)] exp[a(x 1 -1/2)] + exp[-a(x 1 -1/2)] exp[a(x 2 -1/2)] exp[a(x 2 -1/2)] + exp[-a(x 2 -1/2)] ,
x = (x 1 , x 2 ) ∈ X . For a large enough (we take a = 30 in the numerical experiments below), Z(x) is approximately equal to Z 2 (x) for x 1 and x 2 larger than 1/2 and is approximately equal to Z 1 (x) otherwise. The experimental design ξ n = (x 1 , . . . , x n ) is a random Latin hypercube with n = 20 points in X , see Figure 1. Two predictors are compared. The first one is the Empirical BLUP (EBLUP). A stationary GP model is assumed, with unknown constant mean β, unknown variance σ 2 and covariance function K γ 1 ,θ (•) with unknown θ. The parameters β, σ and θ are estimated by Restricted ML, see Section 5.1; these estimated values are then plugged in the ordinary Kriging predictor (the BLUP) -it corresponds to (21) in the Appendix, with g(x) ≡ 1, = 1 and no conditioning on t. We denote by ηMV,n (t) this prediction at t.

The second one is the BLK predictor ηBLK,n (t) of Section 2.2, with L = 4 and covariance functions given by (1) with k 0, (x-x ) = K 3/2,θ 0, ( x-x ) and k 1, (x-x ) = K 3/2,θ 1, ( x-x ), = 1, . . . , 4. We suppose that g(x) ≡ 1 and that β has a uniform prior, see Section 5.1; we take γ 0 = 2, σ 0 = 1 which corresponds to a very vague prior on σ 2 .

We construct predictions on a regular grid of 21 × 21 = 441 points s i in X and consider different values for the parameters of the process ( 16):

θ 1 = 1, β1 = 0, σ 1 = 1, γ 1 = 3/2 or 5/2, γ 2 = 3/2, θ 2 = 7, σ 2 = 2,
β2 = 0 or 10. Note that the EBLUP assumes that the covariance is K γ 1 ,θ (•). BLK uses θ 0, = 1, 5, 10, 20 and either θ 1, = 0 (i.e., it assumes stationarity) or θ 1, = θ * > 0 for all = 1, . . . , 4 (non-stationarity). The value θ * is chosen according to the following rule. Assuming that the design ξ n is space-filling, we wish to ensure that, for any x ∈ X , k 1, (x -x i ) is large enough (say, larger than 20%) for a significant set of design points x i ∈ ξ n (say, 10d/4 points). For the covariance Matérn 3/2, we obtain that the hypercube with side length τ * = 2.9943/θ * should contain 10d/4 points, which yields τ * = min{1, [10d/(4n)] 1/d }.

For n = 10d and d = 2, this gives τ * = 0.5 and θ * 5.9886. The different configurations considered are indicated in Table 1. In columns A-C the process Z(x) is stationary: in A and B the EBLUP uses the correct covariance function, BLK does not assume stationarity in A but does in B and C; the EBLUP has a wrong covariance function in C. The process Z(x) is nonstationary in D and E, with Z 1 (x) and Z 2 (x) having different covariances in D and also different means in E. Figure 2 presents realizations of the process Z(x) for configurations corresponding to columns D and E of Table 1. For each choice of covariance structure for the process Z(•) we repeated 1,000 independent simulations of the 20-dimensional vector of observations z n . To compare the performances of the two predictors without simulating realizations of Z(s i ) for the 21 × 21 grid points, one may notice that, for any predictor ẑn (x 0 ) which is a function of z n , we have

E 2 [ẑ n (x 0 )] = E{[ẑ n (x 0 ) -Z(x 0 )] 2 |z n } = ∆ 2 n [ẑ n (x 0 )] + V n (x 0 ) , with ∆ 2 n [ẑ n (x 0 )] = [ẑ n (x 0 ) -E{Z(x 0 )|z n }] 2 and V n (x 0 ) = E{[Z(x 0 ) -E{Z(x 0 )|z n }] 2 |z n }
, where E{Z(x 0 )|z n } and V n (x 0 ) can easily be calculated using the characteristics of Z(•), see Remark 1. The values of the integrated squared errors

{IE 2 [ẑ n ]} k = 1 441 441 i=1 E 2 [ẑ n (s i |z (k) n )]
k are then calculated for each vector of simulated data z

(k)
n , for both predictors ηMV,n and ηBLK,n . The empirical means

IE 2 [ẑ n ] = 1 1, 000 1,000 k=1 {IE 2 [ẑ n ]} k
of the integrated squared errors are indicated in Table 1, together with the results of pairedcomparisons tests T pc for the differences

{IE 2 [η M V,n ]} k -{IE 2 [η BLK,n ]
} k and their corresponding p-values (see, e.g., Kanji (1993, p. 30)). A box-plot of these differences is presented in Figure 3, indicating that BLK produces more stable predictions than the EBLUP. We also indicate in

)-IE 2 (η BLK,n ) Figure 3: Box-plot of the differences {IE 2 [η M V,n ]} k -{IE 2 [η BLK,n ]} k for the 5 cases considered in Table 1.
one situation only among those considered (column A of Table 1): the process is stationary, the EBLUP has the correct covariance function and BLK does not assume stationarity. For the same random process, the situation is reversed when BLK makes use of the stationarity assumption (column B); the superiority of BLK over the EBLUP increases when ηMV,n wrongly assumes a Matérn 5/2 covariance (column C). BLK was always found to be superior to the EBLUP when Z(•) is non-stationary, as illustrated by columns D and E, where Z 1 (•) and Z 2 (•) differ by their regularity (column D) and also by their mean (column E). Note that predictions are much more precise in column D than in the others, due to the fact that that the process is quite smooth (γ 1 = 5/2) and has a rather large correlation length (θ 1 = 1) on a big part of X ; see in particular the small values of M edE 2 [η M V,n ] and M edE 2 [η BLK,n ] in columns D and E.

We also indicate in the table the values of the empirical errors predicted by the two modeling approaches. For the EBLUP, the squared prediction error at x 0 is given by σ2 n (x 0 ) ρ 2 n (x 0 ), see ( 19) and ( 23). For BLK with θ 1, > 0, observations that are far away from x 0 have negligible influence on the prediction of Z(x 0 ). We thus construct an "equivalent number of observations" n (x 0 ), given by n (x 0 ) = n i=1 k 1, (x 0 -x i ). This value is substituted for n in ( 18) and ( 19)

for the evaluation of the posterior squared prediction error (13) (but not for the evaluation of the marginal likelihood L (z n | , t)). We then compute the empirical means (over the 1,000 simulations and 441 grid points) of these posterior squared prediction errors, which we denote by IE 2 [η M V,n ] and IE 2 [η BLK,n ], to be compared respectively with

IE 2 [η M V,n ] and IE 2 [η BLK,n ].
It is well known that using the plug-in mean-squared prediction error of the BLUP underestimates the true prediction error, the reason being that the uncertainty due to the estimation of the covariance parameters is not accounted for, see Stein (1999, Section 6.8). The table corroborates this result. On the other hand, the error predicted by BLK slightly overestimates (columns A and B,C) or slightly underestimates (column D) the true empirical error, with an exception for column E where the abrupt change in the mean β yields large and hardly predictable errors (see the values of M axE 2 [η BLK,n ]).

Example 2: expected improvement and deceptive function. Kriging prediction can be used for the global optimization of a function f ( [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF]; [START_REF] Mockus | Bayesian Approach to Global Optimization, Theory and Applications[END_REF], a method which has been popularized under the name of Expected Improvement (EI), see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. The function is considered as the realization of a GP, whose characteristics (in particular the parameters θ of the chosen covariance function) are estimated from observations that correspond to evaluations of the function, generally by ML. However, when the estimated parameters are plugged into the Kriging predictor and its associated mean-squared prediction error, the performance of the method may be rather disappointing since evaluation results may not contain enough information to estimate the covariance parameters in a satisfactory manner. This may wrongly provide the sensation that the function is extremely flat in some areas, that will thus not be explored, or on the opposite extremely wiggly, so that all X would seem to deserve a close exploration. This phenomenon is well described in [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF] through the concept of deceptive function, an example of which is presented below. Consider the function in (Benassi et al., 2011, Section 5.1), f (x) = x [sin(10x + 1) + 0.1 sin(15x)], x ∈ X = [-1, 1], plotted in solid line in Figure 4. When the function is evaluated at the design ξ n = (-0.43, -0.11, 0.515, 0.85) all observations are nearly zero, see the stars in Figure 4. We assume that f (•) can be represented by a GP with unknown mean β and Matérn 3/2 covariance function, see ( 15). The EBLUP (blue solid line) and associated (approx-imate) 95% credible intervals (dashed lines), corresponding to ηMV,n (x 0 ) ± 3.182 σn (x 0 )ρ n (x 0 ), with 3.182 the critical value of the student t-distribution with n -p = 3 degrees of freedom, see Santner et al. (2003, p. 95), are plotted in Figure 4-left. The most promising region in terms of maximization of f (•) is around zero and several additional evaluations of f (•) are required before the method starts exploring the neighborhood of the maximizer of f (•), at x * -0.9052, see [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF].

•) : x ∈ X ⊂ R d -→ f (x) ∈ R, see
Figure 4-right presents the BLK predictor and approximate 95% credible intervals given by ηBLK,n (x 0 ) ± 2.571 [P SP E(η BLK,n (x 0 ))] 1/2 (with 2.571 the critical value of the student tdistribution with n -p + ν 0 = 5 degrees of freedom), under the same setting as in Example 1 and θ 0, = 1, 5, 10, 20, θ 1, = 5 for = 1, . . . , L = 4. The large uncertainty on the behavior of f (•) in the left-hand side of the domain is an incitation to put observations there. The EI algorithm can be expected to perform much more efficiently when based on BLK, even with a small L, than when based on the EBLUP.

Example 3: behaviour of BLK far from design points. This example illustrates the behavior of the BLK predictor ηBLK,n (t) when t is far enough from design points so that the influence of the closest point x i dominates all others through k 1, (t -x i ), see (1).

We consider a process Z(x) with zero mean and covariance K 3/2,2 (•), see ( 15), which is observed at the 7 design points (-1, -0.9, -0.8, -0.5, -0.2, 0, 1), indicated by stars. Figure 5 shows the predictions obtained by ordinary Kriging (with covariance K 3/2,15 (•)) and BLK (with L = 4, k 0, (δ) = K 3/2,θ 0, (|δ|) and θ 0, = 1, 5, 10, 20, k 1, (δ) = K 3/2,5 (|δ|) for all ), both assuming that g(x) ≡ 1. The prediction by limit-Kriging [START_REF] Joseph | Limit kriging[END_REF], with covariance K 3/2,15 (•), is also presented on the figure. The correlation length θ being underestimated (1/15 instead of 1/2), the ordinary-Kriging prediction ηn (t) is close to the estimated mean βn given by (20) for t far enough from design points, see the right part of the figure. On the other hand, in the same region the BLK predictor is mainly influenced by the closest design point x i , so that ηn (t) is close to Z(x i ), a reasonable behavior that resembles that of limit-Kriging with underestimated θ, see Joseph (2006, Section 3).

Example 4: prediction of an oceanographic field, using outputs of a formal (numerical) model. The data used in this study were made available through a collaboration with the institute MUMM, a department of the Royal Belgian Institute of Natural Sciences. It consists of snapshots of the output of the biogeochemical oceanographic model MIRO&CO [START_REF] Lacroix | Validation of the 3D biogeochemical model MIRO&CO with field nutrient and phytoplankton data and MERISderived surface chlorophyll α images[END_REF], run to simulate the evolution of inorganic and organic carbon and nutrients, phytoplankton, bacteria and zooplankton with realistic forcing conditions. The model covers the entire water column of the Southern Bight of the North Sea. In the study presented here we concentrate on a surface grid G of 49 × 21 points, 879 of which corresponding to sea surface and form our design space X . The objective is to assess the possibility of predicting chlorophyll concentration over X from observations at a small number of sites. The simulation model is used as a substitute for the real phenomenon; it provides pseudo-observations at n design points ξ n = {x 1 , . . . , x n } in X (n = 25) and allows the computation of empirical prediction squared errors.

Experimental design. Figure 6-Left presents the model response f (x), x ∈ X . It is manifest that variability is stronger along the French coast, so that obtaining precise predictions there would require a denser concentration of observation sites than in other areas where the response is more flat. However, in a realistic situation the true response is not available, and this information cannot be used to choose ξ n . The design ξ n is thus constructed as follows. First, G is renormalized to [0, 1] 2 and X is renormalized accordingly. Then, we generate a low-discrepancy (Sobol') sequence in G ; the design points x 1 , . . . , x 10 are given by the first 10 points of the sequence that fall in X . Finally, the next 15 points are generated sequentially, according to

x k+1 = arg max x∈X ρ k (x) , (17)
with ρ k (x) the universal Kriging variance (the mean-squared error of the BLUP) for a process with unknown constant mean β and covariance function K 3/2,θ (∆(x, x )), see ( 15), and observations f (x 1 ), . . . , f (x k ). We take θ = 10 to ensure that ξ n will be well dispersed over X . In order to take the non-convexity of the design region into account, the "distance" ∆(x, x ) corresponds to the shortest maritime route between x and x , computed by Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]. The 25 design points x i of ξ n are indicated by stars in Figure 6-Left. Figure 6-Right presents the distance ∆(x, x 21 ) from x ∈ X to the 21st design point x 21 (with indices (15, 47), corresponding to coordinates (0.2593, 0.9583) in the renormalized space), showing a neat difference with respect to the Euclidean distance x -x 21 . Notice that in order to compute predictions and prediction errors over X we only need to compute distances ∆(x i , x) from the design points x i to the x ∈ X (and not all pairwise distances ∆(x, x ), (x, x ) ∈ X 2 ).

Comparison between BLUP and BLK For this same design ξ n , we compare the predictive performance of the EBLUP and BLK. The EBLUP uses the covariance function K 3/2, θn (∆(x, x )), with θn 19.2058 estimated by ML from the observations f (x 1 ), . . . , f (x n ); BLK uses K 3/2,θ (∆(x, x )) both for k 0, and k 1, , with L = 4 and θ 0, = 5, 10, 20, 30 in k 0, , θ 1, = θ * > 0 for all = 1, . . . , 4 in k 1, . The choice of θ * should make a compromise between stepping the influence of distant points down (which means taking θ * large enough to obtain local predictions) and maintaining a reasonable correlation with sufficient design points (which means taking θ * not too large). Let φ mM (ξ n ) = max x∈X min x i ∈ξn x -x i denote the value of the minimax distance criterion for ξ n , and denote φ kN N (ξ n ) = max x i ∈ξn x i -x j * (i) with x j * (i) the k-th nearest-neighbor of x i in ξ n . For all x ∈ X , we can then guarantee that

x -x i ≤ δ k (ξ n ) = φ mM (ξ n ) + φ kN N (ξ n )
for k + 1 points x i of ξ n . We take k = 19 and θ * = 2.99431/δ k (ξ n ), which ensures that, for each x ∈ X , k 1, ( x -x i ) ≥ 20% for at least 20 design points in ξ n . For the design ξ n plotted in Figure 6-Left, this gives θ * 2.4202, the value we use here for BLK.

Figure 7-Left shows the BLK predictions over X , to be compared with the true responses on Figure 6-Left. Figure 8 presents the difference in squared prediction errors between the EBLUP and BLK, E 2 [η M V,n (x)] -E 2 [η BLK,n (x)] for x ∈ X . The accuracies of the BLK predictor and the BLUP are similar in most of the domain X , but BLK is more accurate in a large portion of the French coast, precisely where good predictions are difficult to obtain. The mean, median and maximum values of E 2 [η M V,n (x)] over X are respectively 0.0144, 0.0035 and 0.6447; these values equal 0.0122, 0.0035 and 0.5070 for BLK; a paired-comparisons test for these squared errors gives approximately 5.52, with an associated p-value 1.7 10 -8 indicating that BLK is significantly more accurate than the EBLUP. Limit-Kriging [START_REF] Joseph | Limit kriging[END_REF] with θ estimated by ML yields mean, median and maximum values of squared prediction errors respectively equal to 0.0140, 0.0035, 0.6527, and thus performs similarly to the EBLUP. The behavior of the EBLUP is marginally improved when the ML estimation of θ in K 3/2,θ (∆(x, x )) is replaced by leave-oneout cross validation, see [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]: the estimated θ is then θn 17.64, the mean, median and maximum values of squared prediction errors over X become respectively 0.0142, 0.0036 and 0.6312; the paired-comparisons test with BLK gives 5.82, with a p-value 3.0 10 -9 .

[η M V,n (x)] -E 2 [η BLK,n (x)].
Choice of θ * . The influence of the choice of θ * on the performance of BLK (with L = 4 and θ 0, = 5, 10, 10, 30 in k 0, ) is illustrated in Figure 9, where the solid line (respectively dashed line) gives the mean (respectively 1/40 × the maximum) value of E 2 [η BLK,n (x)] over X as a function of θ * varying between 0 (stationary model) and 10 (strong non-stationarity). The choice θ * = 2.99431/δ 19 (ξ n ) 2.4202 is not optimal but seems reasonable. Imposing that each x ∈ X has only at least 10 neighboring design points x i such that k 1, ( x -x i ) ≥ 20% would have been a better choice, since 2.99431/δ 9 (ξ n ) 4.2531, closer to the optimum θ * in Figure 9. Note that the errors E 2 [η BLK,n (x)] are normally not available, and therefore cannot be used to select θ * . This indicates, however, that the choice of k in the construction θ * = 2.99431/δ k (ξ n ) is not critical.

Influence of L. Increasing L does not necessarily improve the performance of BLK. For instance, taking θ 0, = 5, 6, 7, . . . , 30 in k 0, (L = 26) and θ 1, = 2.4202 for all = 1, . . . , 26 in k 1, , we obtain 0.0127, 0.0035, 0.5388, respectively for the mean, median and maximum values of E 2 [η BLK,n (x)] over X . Figure 10-Left shows the posterior weights w n (x) associated with the correlations lengths 1/θ 1, , = 1, . . . , L = 4, when θ 0, = 5, 10, 20, 30, pointing out areas where the response exhibits strong variability and those where it is fairly smooth. It should be stressed that the calculation of the w n (x) only uses the 25 response values f (x i ) for x i ∈ ξ n . The similarity between the maps of w 3 n (x) and w 4 n (x) is an incitation to try reducing L to 3: when θ 0, = 5, 10, 30 in k 0, (L = 3), we obtain 0.0122, 0.0035 and 0.5040 for the mean, median and maximum values of E 2 [η BLK,n (x)]; i.e., values that are marginally better than those with L = 4.

Influence of the distance function. Finally, the interest of using the distance ∆(x, x ) corresponding to the shortest maritime route between x and x instead of the Euclidean distance x-x (in the renormalized space [0, 1] 2 ) is illustrated for BLK in Figure 10-Right, which shows the differences between corresponding prediction squared errors. The distance ∆(x, x ) yields more accurate predictions along the major part of the French coast; the mean, median and maximum values of E 2 [η BLK,n (x)] over X are respectively 0.0129, 0.0030 and 0.5257 when using Euclidean distance to compute correlations. The decrease of performance compared to the situation where ∆(x, x ) is used is not caused to the fact that the sequential construction of design points x 11 to x 25 through ( 17) is based on ∆(x, x ): replacing the design ξ n by ξ n where x 11 to x 25 are constructed via (17) with Euclidean distance based covariances, we obtain mean, median and maximum squared prediction errors equal to 0.0142, 0.0034, 0.5175 when BLK uses Euclidean distance.

Conclusions and further developments

We have presented a Bayesian local Kriging approach for the interpolation or prediction of random fields. The method uses localized covariance functions which allow us to account for non-stationarity. A finite set of L candidate covariance functions is used for each prediction point, that receive equal prior weights. Using a common hierarchical prior for the trend and variance of the process, posterior weights can easily be calculated to obtain posterior means and prediction squared errors. Numerical simulations indicate that, on average, the method performs slightly worse than universal Kriging with plug-in maximum-likelihood estimates for the covariance parameters when the true characteristics of the process satisfy the assumptions (stationarity, correct parametric trend, correct covariance functions), but performs significantly better when these assumptions are violated, even if the number L of concurrent covariance functions is very small (L = 4 in the examples considered). Also, BLK, which does not use any numerical optimization, is much faster and numerically stable than universal Kriging with plug-in estimates, which requires the estimation of covariance parameters.

To summarize, our feeling is that it seems illusory in many applications to try estimate covariance parameters from a few observations only, especially with a covariance structure not necessarily well-adapted to the variability of the modeled phenomenon. Using a small number of candidate processes able to reproduce a reasonable range of possible behaviors may be preferable: it is simpler to implement, numerically more stable, and seems to often yield better predictions.

Although these results are encouraging, further numerical experimentations (in particular for higher dimensional processes and different types and sizes of experimental designs) are needed to confirm these preliminary observations. We have restricted our attention to the situation where the regressor g(•) in the parametric trend g (x)β was identical for all L models (and, moreover, we only considered the case g(x) ≡ 1 in all examples). The same approach could be used when different trends g (•) are associated with different covariances C |t (•, •), possibly allowing for a better consideration of uncertainty in the process trend.

The choice of the particular form of the kernels k 0, (•) and k 1, (•) seems to be less crucial than that of the correlation length for k 1, (•), which depends on the assumed amount of nonstationarity. Using different correlations lengths for some of the L concurrent covariances is a possible option to investigate. Another one is to simply ensure that k 1, (x -x i ) be large enough for all points in x in X and enough design points x i , with the motivation that the more dense the design ξ n in X , the more local the models can be and the stronger the non-stationarity that BLK can take into account. Two proposals have been made in this direction, see Examples 1 and 4. Further investigations are required to validate them, in particular concerning their asymptotic behavior, when the number n of design points tends to infinity and ξ n is space-filling.

Finally, as mentioned in Section 2.3, designing experiments adapted to BLK requires the evaluation (approximation) of the second term in ( 14). This rather challenging problem is under current investigation.

5 Appendix: Universal BLK in presence of a parametric trend Here we consider that, seen from t, f (•) is the sample path of a stochastic process Z(•) such that Z(•)|s(t), β, σ 2 , t ∼ GP(g (•)β, σ 2 C s(t)|t (•, •)), with Prob{s(t) = } = w 0 , = 1, . . . , L, β|σ 2 ∼ ψ 0 (β|σ 2 ), σ 2 ∼ ϕ 0 (σ 2 ), see (8), and g(•) a known p-dimensional vector of functions defined on X . The likelihood of z n = (f (x 1 ), . . . , f (x n )) for the model GP(g (•)β, σ 2 C s(t)|t (•, •)) is

L (z n |β, σ 2 , , t) = 1 σ n (2π) n/2 det 1/2 K n ( |t) exp - 1 2σ 2 (z n -G n β) K -1 n ( |t)(z n -G n β) ,

Figure 1 :

 1 Figure 1: Experimental design ξ n in Example 1.

Figure 2 :

 2 Figure 2: Realizations of a non-stationary process Z(x) with w(x) = w n-s (x), for configurations D (left)and E (right) of Table1.

Figure 4 :

 4 Figure 4: Left: a deceptive function (solid line), with the EBLUP (dashed line) and 95% credible intervals (dotted lines). Right: the same function with BLK (dashed line) and 95% credible intervals (dotted lines).

Figure 5 :

 5 Figure 5: Predictions for BLK (dashed line with triangles), ordinary Kriging (dash-dotted line with squares) and limit Kriging (dotted line with circles) with underestimated correlation length; Z(x) is in solid line, observations are indicated by stars.

Figure 6 :

 6 Figure 6: Left: model response over X and design points (stars). Right: shortest maritime distance ∆(x, x 21 ) between x ∈ X and the 21st design point x 21 .

  Figure 7-Right shows the squared prediction errors E 2 [η BLK,n (x)], where for a predictor ẑn (•) we denote E 2 [ẑ n (x)] = [ẑ n (x) -f (x)] 2 . Taking into account that only 25 designs points have been used, predictions are fairly accurate, excepted at some areas along the French coast where the correlation structure strongly departs from the smoother variation in the open sea region.

Figure 7 :

 7 Figure 7: Left: BLK predictions ηBLK,n (x) over X . Right: squared errors E 2 [η BLK,n (x)] for BLK.

Figure 8 :

 8 Figure 8: Difference in squared errors between the EBLUP and BLK, E 2 [η M V,n (x)] -E 2 [η BLK,n (x)].

Figure 9 :

 9 Figure 9: Mean value (solid line) and 1/40 × maximum value (dashed line) of E 2 [η BLK,n (x)] over X as functions of θ * .

Figure 10 :

 10 Figure 10: Left: posterior weights w n (x) for correlations lengths 1/θ 1, . Right: differences between squared errors for BLK with Euclidean distance x -x and BLK with ∆(x, x ).

Table 1 :

 1 Table 1 the values M edE 2 [ẑ n ] and M axE 2 [ẑ n ] of the empirical means, over the 1,000 simulations, of respectively the median and maximum values of E 2 [ẑ n (s i |z The values of the paired-comparisons tests T P C and associated p-values in Table 1 indicate that conclusions about the best predictor (in terms of integrated squared errors) between ηMV,n and ηBLK,n are highly significant. The EBLUP, that is, ordinary Kriging with ML estimation of the process parameters, appears to yield more precise predictions (on average) than BLK in Simulation results in Example 1. Best performances (smallest values) among ηMV,n and ηBLK,n are indicated in bold face.

	(k) n )]
	8

k on the 441 grid points.

* This work was partly supported by the ANR project 2011-IS01-001-01 DESIRE (DESIgns for spatial Random fiElds), joint with the Satistics Department of the JKU Universität, Linz, Austria.

where the i-th row of the n×p matrix G n equals g (x i ), i = 1, . . . , n. The weights w l 0 are updated according to (3), with now L (z n | , t) = R p ∞ 0 L (z n |β, σ 2 , , t) ψ 0 (β|σ 2 ) ϕ 0 (σ 2 ) dβ dσ 2 .

Uniform prior for β

Suppose first that β has a uniform (improper) prior on R p . We obtain

is the Restricted Maximum-Likelihood (REML) estimator of σ 2 , and

is the ML estimator of β, given z n , and t, see, e.g., Santner et al. (2003, p. 67, 95). Moreover, given z n , and t, σ 2 has the inverse chi-square distribution Santner et al. (2003, p. 95)

and βn ( |t) given by (20). Therefore, ηn (x 0 |t) = c n (x 0 | , t)z n and ηn (x 0 |t) is still given by ( 11), but with

with I n the n-dimensional identity matrix. Also,

or equivalently,

Therefore, the posterior squared prediction error is still given by ( 13), but with ν n = ν 0 + n -p, σ 2 n| ,t given by ( 18), ρ 2 n (x 0 | , t) by ( 23) and c n (x 0 | , t) by ( 22) in Ω n (x 0 |t). Similarly to Sect. 2.3, the preposterior variance at t is (for ν 0 > 2) is given by ( 14).

Normal prior for β

Suppose now that ψ 0 (β|σ 2 ) is the density of the p-dimensional normal distribution with mean β 0 and variance-covariance matrix σ

so that, similarly to Sect. 2.3,

where ν n = ν 0 + n and

see ( 24), with now

or equivalently,