Egor Sattarov

Alexander Gepperth

Sergio A Rodriguez

Roger Reynaud

Calibration-free match finding between vision and LIDAR*

We present a learning approach that allows to detect correspondences between visual and LIDAR measurements. In contrast to approaches that rely on calibration, we propose a learning approach that will create an implicit calibration model from training data. Our model can provide three functions: first of all, it can convert a measurement in one sensor into the coordinate system of the other, or into a distribution of probable measurements in case the transformation is not unique. Secondly, using a correspondence observation that we define, the model is able to decide if two visual/LIDAR measurements are likely to come from the same object. This is of profound importance for applications such as object detection or tracking where contributions from several sensors need to be combined. We demonstrate the feasibility of our approach by training and evaluating our system on tracklets in the KITTI database as well as on a small set of real-world scenes containing pedestrians, in which our method finds correspondences between the results of real visual and LIDAR-based detection algorithms.

I. INTRODUCTION

A. Context of this work

This article is in the context of multisensory information processing, in particular vision and LIDAR. As these sensors take their measurements independently, it is a priori not clear whether two measurements originate from the same object (or more generally: from the same physical position). To find these correspondences, standard algorithms like object detection and tracking (i.e. DATMO) usually make use of a calibration procedure which allows to transform measurements of one sensor into the reference frame of the other. Such transformations are often quite sensitive [START_REF] Rodriguez | Influence of intrinsic parameters over extrinsic calibration between a multi-layer lidar and a camera[END_REF] to the used measurement models (e.g., pinhole model for camera) and calibration parameters. Moreover, because of the very nature of the measured quantities, sometimes a one-to-one transformation does not even exist. This is for example the case when transforming 2D image points into a 3D coordinate system of a LIDAR device.

B. Proposed approach

While calibration approaches are often quite precise, the calibration procedure itself is complex and error-prone and requires considerable expertise. Furthermore, a calibration procedure intrinsically depends on the common data representation (e.g. calibration pattern, features), and needs *This work was supported by a PhD grant from the Digiteo consortium 1 Egor Sattarov, Sergio A. Rodriguez F. and Roger Reynaud are with Université Paris-Sud and the Institute d'Electronique Fondamentale UMR 8622 CNRS, Orsay, France {egor.sattarov, sergio.rodriguez and roger.reynaud}@u-psud.fr to be re-designed every time a change is made. On the other hand, it is often rather easy and cheap to obtain a large number of sample measurements from both sensors. Assuming the existence of such a sample database, we propose a simple method to extract an implicit calibration model between vision and LIDAR sensors. We pursue a data-driven approach where the statistics of each sensor are used to optimally project both measurements (i.e. objectlevel) onto a standardized representation format to which generic probabilistic methods can be applied. In this way, our approach is completely independent of the intrinsic characteristics of the measurements, and in particular of their dimensionality (i.e. n-d observations), leading to a strong reduction of design and re-design effort for the conception of multi-modal processing system in vehicles.

C. Related work

State-of-the-art approaches for IV such as [START_REF] Cho | A multisensor fusion system for moving object detection and tracking in urban driving environments[END_REF], [START_REF] Rodriguez | Multimodal object detection and localization for high integrity driving assistance[END_REF] rely on the explicit need of a common frame where all sensors observations can be referenced (i.e. data alignment). This assumption greatly simplifies the association problem of multiple data sources (e.g. LIDAR, Radar, vision). However, in practice, a calibration procedure is required in order to precisely determine all sensors rigid-body transformations (i.e. extrinsic parameters) into the reference frame and their uncertainties.

Recent works on the 3D sensor calibration have considerably simplified the procedure for determining the relative position of sensors using a set of natural features [START_REF] Scaramuzza | Extrinsic self calibration of a camera and a 3d laser range finder from natural scenes[END_REF] or using a single observation of a set of calibration patterns ( covering different distances and orientations of the multi sensors field of view) [START_REF] Geiger | Automatic calibration of range and camera sensors using a single shot[END_REF].

Automatic calibration approaches can also infer the extrinsic parameters by the means of an optimization framework which registers sensors data in a common space (typically 2D/3D Cartesian space). Recently in [START_REF] Pandey | Automatic extrinsic calibration of vision and lidar by maximizing mutual information[END_REF], [START_REF] Levinson | Automatic online calibration of cameras and lasers[END_REF] and [START_REF] Napier | Cross-calibration of push-broom 2d lidars and cameras in natural scenes[END_REF], online strategies were proposed to achieve data registration between a vision system and a ranging sensor by optimizing the extrinsics following a mutual information criterion of the sensing sources.

As an alternative to the classical approaches , the presented work is intended to perform multi-sensor data alignment through a probabilistic learning based framework. This approach not only provides a data alignment solution but also models the probability accorded to the observation transfer. Moreover, this method can provide an integrity measure of the data alignment using extrinsic parameters in a crossvalidation scheme.

D. Contributions and novelty

This article presents a new way of detecting multimodal correspondences for the important vision/LIDAR sensor combination that is becoming a standard in the intelligent vehicle domain. A main contribution of the used learning approach is that the "calibration" procedure is much simpler and can in fact be handled by a non-expert regardless of the precise type of measurements that are conducted. Furthermore, we show that the resulting data alignment is very computationally efficient and sufficiently accurate for most applications. Performing all experiments using the publicly available KITTI benchmark database adds significant credibility to our results.

II. METHODS

A. Architecture overview

The complete model is composed of several components, as visualized in Fig. 2: • LIDAR and vision sensors • means to measure interesting quantities in both • Self-organized Maps (SOM) for vision and for LIDAR, which learn to represent the inputs coming from the respective (synchronous) measurements • an algorithm for learning a correspondence model between SOMs • a module for deciding when two measurements correspond, based on the SOMs and the learned correspondence model Within the scope of this article, we will use both actual measurements proposed by Honda experimental vehicle in form of object positions, and annotated tracklets from the KITTI database [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF] as ideal data without noise. They are used separately. As we wish not to complicate the clean and simple algorithm we propose by details of unimodal processing in each modality.

B. Model training 1) Learning sensor statistics with self-organizing maps:

The self-organizing map algorithm, while originally proposed as a model cortical information processing, is a generative machine learning algorithm that aims to approximate the distribution of high-dimensional data, and to represent it in a topology-preserving way on a two-dimensional manifold. It is in fact quite related to K-Means [START_REF] Jain | Algorithms for Clustering Data[END_REF] except that the preservation of topology makes it interesting for incremental learning scenarios.

SOM defines a fixed N × N grid of nodes ("neurons") n i , each of which is associated with a so-called prototype vector p i . For a given input x, each node gets assigned an activity z i based on the distance of its prototype to the input:

z i = d( x, p i ) d( a, b) = ( a -b) 2 (1)
As a distance measure, the euclidean distance is often used, and so shall we. In most cases, the calculation of activity is followed by a learning step where the prototypes are adapted to better fit the current input:

i * = argmin i z i p i (t + 1) = p i + (t)G (i * , i, σ(t)) ( p i * -p i ) (2)
G(i, j, σ = exp (-d 2 (i,j) 2σ 2
)) is a Gaussian with standard deviation σ which is based on the euclidean distance between node i and node j on the two-dimensional grid of nodes. For faster convergence, the algorithm demands to gradually lower the learning rate (t) and their neighbourhood radius σ(t) from initially large values 0 , σ 0 until the minimal values ∞ , σ ∞ are reached.

2) Learning of conditional distributions between sensors: Supposing the SOMs are trained using the algorithm described in Sec. II-B.1, correspondences between visual and LIDAR SOMs are detected using a simple probabilistic counting approach. Assuming that two sets of weights w L ij , w V ij exist between nodes i, j in visual and LIDAR SOMs, both are updated as follows for each simultaneously presented pair of visual and LIDAR measurements x V , x L :

zX i = 1 if i = argmin k z X k 0 else z X i =        1 if i = argmin k z X k 0.5 if i is neighbour to argmin k z X k 0 else w X ij = w X ij + zX i z X j (3)
where we have used a shorthand notation X = L, V ( X denoting the other modality, i.e., L if X = V and V otherwise). After a sufficient amount of samples has been processed, we normalize the weight matrices in order to obtain normalized probabilities:

Σ X i = j w X ij w X ij → w X ij Σ X i (4)
It has to be noted that the visual and LIDAR measurements do not need to come from the same objects. Indeed, if this were the case, it would mean that we already know the correspondences we wish to identify. When working on a benchmark database like KITTI, this is the case but when training the system on recorded data not containing any annotations, evidently the correct correspondences are unknown except when there is always just a single object in sight. Therefore, the strategy we adopt is to present all combinations of visual and LIDAR measurements taken at a certain point in time (e.g. a single, synchronized image and LIDAR recoding, both for real sensors and in the case for KITTI) when learning conditional probability distribution between sensors. This assumes there is a sufficient amount of training data, because the "correct" correspondences will appear together far more often than random incorrect ones. As we supposed that SOMs are already converged, we disable SOM learning during the whole phase of learning conditional distributions by setting (t) ≡ 0 for both SOMs.

3) Overall training procedure: The overall training procedure is given in Alg. 1. It consists of a SOM training step and a step that determines conditional probabilities between the SOM representations of both measurements.

C. Unimodal detection of correspondences

After training is completed, the model can be used for detecting whether a given combination of visual and LI-DAR measurements is likely caused by the same object. To this end, we develop a criterion that depends on a single parameter, the probability threshold θ. Assuming that each measurement has generated activities z X i in both SOMs, the criterion first computes a single binary measure c X = {0, 1} for each conditional probability matrix w X ij , using the shorthand notation X = L, V for a certain modality, and Draw a random LIDAR measurement m x L im from i

i * = argmin i z X i (5) j * = argmin j z X j (6) P X = j|w X i * j > θ (7) c X = 1 if j * ∈ P X 0 else (8)
The two quantities c X express whether a best-matching unit (BMU) at position i * in X can predict the best-matching unit at index j * in the other modality X based on the learned conditional probabilities. Given a best-matching unit in X, θ is used for selecting a set of nodes P X with conditional probabilities that exceed θ. If the BMU of X is an element of the selected set, we conclude that there is a match and set c X = 1. Thus, the threshold θ governs the strictness of the matching: if it is high, only a small (or empty) set of nodes P X will be selected and the probability of match diminishes. On the other hand, if θ is low, the probability of match increases, up to the point where there will always be a match at θ = 0. As it is often not necessary to detect all correspondences correctly but rather to exclude unlikely combinations, a more relaxed value of θ can help avoid missed correspondences while still being able of reducing the combinatorial space of correspondences.

D. Fused correspondence detection

Apart from the unidirectional mutual sensor activity predictions, one can also use a cross-verified decision of improved quality. For that, the criterion of acceptance in eqn.(8) changes to:

w X i * j * × w X j * i * > θ (9)
w X i * j * + w X j * i * > θ (10) k w X i * k z X k × k w X j * k z X k > θ (11)
where z X i is again the activity at node i in SOM X (which can be LIDAR or vision, whereas X represents the other modality), and the indices i * , j * are the indices of the BMU's in both sensor's SOMs. The last eqn. [START_REF] Heskes | Energy functions for self-organizing maps[END_REF] takes into account not only the BMU of each SOM, but also its neighbouring nodes plus their associated, learned conditional probabilities.

E. Training and evaluation data

For the training stage and the evaluation of the proposed methods, two datasets are used:

• Dataset A is composed of annotated tracklets from the public KITTI benchmark database [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF] (see also Fig. 5). • Dataset B is composed of real detections captured from dash camera and four-layer lidar on-board an experimental platform. Visual pedestrian detections are obtained with the 'daimer' detector provided with the OpenCV vision library, and LIDAR ones are from connectivity based clustering. From Dataset A: we employ the center positions of objects in 2D image coordinates as well as corresponding 3D laser coordinates as measured by a Velodyne laser scanner (i.e. tracklets). As the height-over-ground of a tracklet's center is often irrelevant for safety applications, we take a birds-eye perspective and just consider two of the three 3D coordinates, excluding height-over-ground. Due to the synchronized nature of visual and LIDAR recordings in the Dataset A, each tracklet can be assigned a unique visual image and therefore a corresponding LIDAR sweep. We use all types of objects provided by the database A, making the total number of considered tracklets 23497. For training the model, we use 70% of this data, performing a random split of available tracklets into train and test databases.

From Dataset B: we employ the center positions of objects in 2D image coordinates as well as corresponding 3D laser point cluster center positions. Vision-based objects and lidarbased detection were manually associated so as to obtain a ground-truth reference. This sequences is composed of pedestrians filmed in 9 short scenarios of about 2 minutes, making a total number of 8613 visual detections, and 5476 lidar detections. Due to the small size of this data base, we performed a cross-validation, that is, for each scenario the SOM are trained with 8 other scenarios and tested with the chosen one.

F. Evaluation

In order to quantify the capacity of the trained model to identify visual/LIDAR correspondences, we use the test database as described in Sec. II-E. In order to prevent the SOMs from adapting during the evaluation phase, we set (t) ≡ 0 for both SOMs. Assuming a trained model (SOMs plus conditional probabilities), we process all images in the test database in a Algorithm 2 Evaluation: Overview over the evaluation procedure. for (l, m) = combinations of measurements do

5: Feed visual SOM with x V il → z V (t) 6:
Feed LIDAR SOM with x V im → z L (t)

7:

Generate bin. measures c X , c X acc. to. Sec. II-C 8:

end for 9: end for 10: Plot precision/recall curves sequential manner. For each image, we present all combinations of visual and LIDAR measurements and compute the scores c L , c V for each combination. A binary decision on the presence of a correspondence is taken according to eqn. [START_REF] Geiger | Automatic calibration of range and camera sensors using a single shot[END_REF]. As this decision depends on a single threshold θ we can re-cast this evaluation in the form of a ROC analysis by varying θ in the interval [0, 1] and measuring the precision/recall rates.

An overview over the complete evaluation procedure is given in Alg. 2.

III. EXPERIMENTS A. Organization of training and evaluation

Model training is performed in two steps: initially, the SOMs are trained independently of one another by drawing random samples from the train database, see Sec. II-E, and adapting each individual SOM according to Sec. II-B.1, with the input vector provided by the unimodal part of the drawn sample. Training parameters are: N = 30, ∞ = 0.01 σ ∞ = 1, 0 = 0.6, σ 0 = N 2 . Neighbourhood radius and learning rate develop according to

σ(t) = max (σ ∞ , σ 0 exp(-λ σ t)) (12) (t) = max ( ∞ , 0 exp(-λ t)) , (13)
with -λ = 0.002 and λ σ = 0.004. SOM training duration is limited to T SOM = 20000 iterations. Subsequently, correspondences are trained according to Sec. II-C for another T corr = 20000 iterations, randomly drawing images from the training database and feeding all possible combinations of visual/LIDAR measurements to the two SOMS as well as updating the two sets of weights w V ij , w L ij based on the resulting SOM activities z X j , X = L, V . Evaluation is conducted according to Sec. II-F by iterating over all images in the test database and measuring precision/recall rates when presenting to the model all possible combinations of visual/LIDAR measurements in each image.

B. Results

For KITTI base we first plot a separate ROC for LIDARvision and vision-LIDAR correspondence detection, given in Fig. 6. As can be expected, the LIDAR-vision-based correspondence detection gives better results, very likely because the vision-LIDAR transformation is one-to-one but not the ROCs for vision-LIDAR (red curve) and LIDAR-vision (green curve) correspondence detection, where laser measurements are augmented by object size. By comparison to Fig. 6, we observe that this irrelevant information is ignored. other way round. We observe as well that performance is acceptable given that no prior knowledge was used at all but it is not an ideal ROC either. In a further experiment, we wish to back the claim made in Sec. I that the proposed method was able to handle arbitrary measurements without requiring explicit models. To this end, we repeat the previous experiment while tracklet width and tracklet height to the laser measurement, bringing up its dimensionality to 4. The ROCs obtained in this way are shown in Fig. 7. We see that the addition of additional information does not impair the ability of our system to detect correspondences. On the other hand, performance is not improved either, because the added information is irrelevant to the transformation to be computed. This experiment therefore shows that our model, due to the learning approach, is able to process very diverse types of measurements, and automatically extracts the information required for finding correspondences. Lastly, we evaluate the three fusion strategies proposed in Sec. II-D, which means that for a pair of visual and LIDAR measurements, there will now be only one decision on correspondence, not two as in previous experiments. The overall performance is shown in Fig. 8 and show that the fused decision outperforms any single unimodal one, boosting the already satisfactory ROCs for fused correspondence detection. The blue "complex" curve represents the cross-verified strategy of eqn. [START_REF] Heskes | Energy functions for self-organizing maps[END_REF]. The cyan "simpleadd" curve corresponds to eqn. [START_REF] Jain | Algorithms for Clustering Data[END_REF], and the violet "simple-mult" one to eqn. [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF]. It is apparent that all fusion methods outperform the unimodal ones (red and green curves). performance even further.

For real sensors the ROCs are calculated using only complex activity predictions [START_REF] Heskes | Energy functions for self-organizing maps[END_REF] as the most effective. The results are seen in Fig. [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF]. One can observe low quality for unidirectional correspondences detections and very high quality for fused one. It can be explained by non-symmetrical detections nature and small number of detected objects per frame.

IV. DISCUSSION AND CONCLUSIONS

We have presented a learning approach to solve the problem of finding visual/LIDAR correspondences and validated its performance on a widely accepted benchmark database. In this section, we will review and justify the components of our model and outline principal conclusions and further research works.

A. Model justification

The hybrid SOM-based architecture we propose here is based on two necessities: first, to have a generic model that will work with any kind of visual/laser measurements. This means that the model must be able to work regardless what is actually measured by each sensor. For a camera, this could be, e.g., pixel position of interest points, but also center position, size and identity if an object detection algorithm is used, or center position, size and speed if tracking is added. By using the self-organizing map architecture, every measurement is down-projected to a 2D image-like representation in a way that is statistically optimal and respects a certain topological constraint that allows to easily visualize and interpret a SOM's activity. For ensuring statistical optimality, we use a variant of the SOM model that has a well-defined energy function [START_REF] Heskes | Energy functions for self-organizing maps[END_REF], which makes it actually very easy to detect measurement outliers that should be ignored.

Secondly, we want a model that will not fail even when the transformation between modalities is not one-to-one in both directions. To this end, we adapted a purely probabilistic approach, on top of the SOM mechanism, that will simply respond by a multi-peaked probability distribution in case there is inherent ambiguity due to non-unique transformations.

B. Discussion of results

As seen in Sec. III, the quality of correspondence finding is very satisfactory given that we did not bring in any specific expert knowledge. In addition, the threshold θ allows us to smoothly change the behavior of the system, from a point where there are few correct correspondences but also few incorrect ones, to a point where there are many correct correspondences but also some incorrect ones. For example, for a multimodal tracking system a higher false positive rate can be acceptable if no correspondences are incorrectly rejected, since tracking can take into account past information and thus correct the occasional incorrect correspondence. Another very encouraging fact is that the quality of correspondence detection can be significantly improved by considering not only both unidirectional correspondences in isolation, but a fusion of both. As a proper fusion should be, it is indeed better-performing than any single contribution to it.

C. Conclusion

We have shown that a learning-based approach can successfully solve the problem of multimodal correspondence detection, in particular between visual and LIDAR sensors. The only prerequisite is a collection of (unlabeled) data which is usually easy to obtain. No expert effort is required at all, and in particular no detailed models of the data acquisition process by the used sensors. The technique is very computationally efficient, and consumes no significant computational load, thus making it suitable for embedded operation. We hope to make this technique even more appealing by better exploiting the structure of conditional probabilities for even better-performing fusion strategies.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Illustration of the multisensory correspondence problem: LIDAR (left) and visual (right) measurements, e.g., provided by independent object detection algorithms, "live" in completely different spaces and are thus very difficult to associate without applying prior knowledge.

Fig. 3 .

 3 Fig. 3. Statistical models of sensory spaces acquired by self-organizing maps (SOM) for visual (left) and LIDAR sensors (right). The points represent the position of SOM prototypes in the space of each sensor. The local density of prototypes is guided by average local density of data points.

Algorithm 1

 1 Model Training: Overview over the two-stage model training procedure consisting of learning distributions with SOMs, and learning multi-sensory conditional probabilities. 1: for t : 1 → T SOM do 2: Draw a random image i from D train 3: Draw random visual measurement l x V il from i 4:Draw a random image i from D train 5:

Fig. 4 .

 4 Fig. 4.Example of conditional probability distributions P X for vision (given a LIDAR node, left) and LIDAR (given a vision node, right). These distributions are used to detect correspondences.

Fig. 5 .

 5 Fig. 5. The KITTI database we use for our experiments is recorded from a moving car equipped with several cameras, a GPS device and most notably a Velodyne LIDAR device.

Fig. 6 .

 6 Fig.6.ROCs for vision-LIDAR (red curve) and LIDAR-vision (green curve) correspondence detection. As can be expected, LIDAR-vision provides slightly better performance as the associated transformation is one-toone.

  Fig.8.ROCs for fused correspondence detection. The blue "complex" curve represents the cross-verified strategy of eqn.[START_REF] Heskes | Energy functions for self-organizing maps[END_REF]. The cyan "simpleadd" curve corresponds to eqn.[START_REF] Jain | Algorithms for Clustering Data[END_REF], and the violet "simple-mult" one to eqn.[START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF]. It is apparent that all fusion methods outperform the unimodal ones (red and green curves).

  : end for 9: Disable learning in SOMs by setting (t) ≡ 0 10: for t : 1 → T corr do

	11:	Draw a random image i from D train
	12:	for (l, m) = all permutations of measurements do
	13: 14: 15:	Feed visual SOM with x V il → z V (t) Feed LIDAR SOM with x L im → z L (t) Update w L ij , w V ij acc. to Sec. II-B.2
	16:	end for
	17:	Normalize w L ij , w V ij acc. to Sec. II-B.2
	18: end for
	"other" for the other one:

6: Update visual SOM with x V il acc. to Sec. II-B.1 7: Update LIDAR SOM with x V im acc. to Sec. II-B.1 8

1 :

 1 Disable learning in both SOMs by setting (t) ≡ 0 2: for i : 1 → (images in D test ) do

	3:	Draw image i from D test
	4:

ROCs for fused correspondence detection in the case of real sensors data. The red "complex" curve represents the cross-verified strategy of eqn. [START_REF] Heskes | Energy functions for self-organizing maps[END_REF]. unimodal ones (red and green curves)


Nadia Morsli
email: nadia.morsli@imag.fr

On asymptotic normality of nonparametric estimate for a stationary pairwise interaction point process

Keywords: Georgii-Nguyen-Zessin formula, nonparametric estimation, pairwise interaction point process, Papangelou conditional intensity, asymptotic normality

We prove the asymptotic normality of nonparametric estimator of pairwise interaction function for a stationary pairwise interaction point process characterized by the Papangelou conditional intensity and observed in a bounded window of a sequence of cubes growing up to R d . Formula for the variance of the resulting estimator can be obtained using Papangelou conditional intensity of the point process. This is a random function satisfying the counterpart of the Georgii-Nguyen-Zessin formula. The proof of the asymptotic normality of the resulting estimator is based on the m n -approximation method in the setting of dependent random fields indexed by Z d where d is a positive integer.

Introduction

A special type of Gibbs processes that are noteworthy both for their abundant use in statistical physics are the pairwise interaction point processes. Several methods of estimation are available for parametric families of pairwise interaction point processes: approximate maximum likelihood ( [START_REF] Ogata | Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure[END_REF], [START_REF] Ogata | Likelihood analysis of spatial point patterns[END_REF]), maximum pseudo-likelihood ( [START_REF] Besag | Point process limits of lattice processes[END_REF]), Monte Carlo likelihood ( [START_REF] Geyer | Likelihood inference for spatial point processes[END_REF]).

In this paper we are concerned with nonparametric statistics for stationary pairwise interaction point processes. They have been introduced by [START_REF] Ripley | Markov point processes[END_REF], [START_REF] Daley | An introduction to the Theory of Point Processes[END_REF] and [START_REF] Georgii | Gibbs measures and Phase Transitions[END_REF]. These provide a large variety of complex patterns starting from simple potential functions (or pairwise interaction function) which are easily interpretable 1 as attractive or repulsive forces acting among points and are of practical importance because of their ability to model a wide variety of spatial point patterns, especially those displaying some degree of spatial regularity. A great deal is understood about pairwise interaction models because they are very natural with respect to the derivation of Papangelou conditional intensity ( [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]).

The goal of the present paper is to establish the asymptotic normality of the nonparametric estimator of the pairwise interaction function for a stationary pairwise interaction point process characterized by the Papangelou conditional intensity. The asymptotic normality in statistics is based on central limit theorems for the sequences of random variables. These classical limit theorems have been extended to the setting of spatial processes. Some results on the central limit theorem and its functional versions are [START_REF] Bass | Law of the iterated logarithm for set-indexed partial sum processes with finite variances[END_REF], [START_REF] Basu | On functional central limit theorem for stationary martingale random fields[END_REF], [START_REF] Dedecker | Central limit theorem for stationary random fields[END_REF], [START_REF] Dedecker | Exponential inequalities and functional central limit theorems for random fields[END_REF] and [START_REF] Machkouri | Kahane-Khintchine inequalities and functional central limit theorem for stationary real random fields[END_REF]. Our proof is based on the m-approximation method (m n → ∞ as n → ∞) in the setting of dependent random fields indexed by Z d where d is a positive integer. We apply a central limit theorem for triangular arrays of stationary m-dependent random fields with unbounded m established by [START_REF] Wang | On the asymptotic normality of kernel density estimators for linear random fields[END_REF]. This result improves a central limit theorem established by [START_REF] Heinrich | Asymptotic behaviour of an empirical nearest-neighbour distance function for stationary Poisson cluster processes[END_REF]. The notion of m-dependence was first introduced by [START_REF] Hoeffding | The central limit theorem for dependent random variables[END_REF]. Central limit theorems for m-dependent random fields has already been considered by many researchers. For example [START_REF] Cheng | Central limit theorems for instantaneous filters of linear random fields[END_REF], [START_REF] Wang | A new condition on invariance principles for stationary random fields[END_REF], [START_REF] Rosén | A note on asymptotic normality of sums of higher-dimensionally indexed random variables[END_REF] [START_REF] Machkouri | A central limit theorem for stationary random fields[END_REF].

The paper is organized as follows. Section 2 sets up the basic tools of the Gibbs point processes in R d . Section 3 contains the main asymptotic results on proposed estimators of stationary pairwise interaction point process. This concerns firstly the asymptotic behaviour of variance and covariance of the kerneltype estimator under certain smoothness conditions. The latter is based on the knowledge of Papangelou conditional intensity and the iterated Georgii-Nguyen-Zessin formula. Secondly, we are able to study asymptotic Gaussianity of proposed estimators and prove it in Section 4. O ⊆ B d be the system of all bounded Borel sets. A point process X in R d is a locally finite random subset of R d , i.e. the number of points N(Λ) = n(X Λ ) of the restriction of X to Λ is a finite random variable whenever Λ is a bounded Borel set of R d (see [START_REF] Daley | An introduction to the Theory of Point Processes[END_REF]). We define the space of locally finite point configurations in

Basic tools

R d as N l f = {x ⊆ R d ; n(x Λ ) < ∞, ∀Λ ∈ B d 0 }, where x Λ = x ∩ Λ. We equip N l f with σ -algebra N l f = σ {{x ∈ N l f : n(x Λ ) = m}, m ∈ N 0 , Λ ∈ B d 0 }, where N 0 = N ∪ {0} = {0, 1, 2, 3, . . .}. The volume of a bounded Borel set Λ of R d is denoted by |Λ|.
For any finite subset I of Z d , we denote |I| the number of elements in I, a, b ≡ {a, a + 1, . . . , b} for a, b ∈ Z and the sum ∑ = signifies summation over distinct pairs. For any Gibbs point processes in a bounded window, the conditional intensity at a location u given the configuration x is related to the probability density f by

λ (u, x) = f (x ∪ {u}) f (x) (for u / ∈ x)
, the ratio of the probability densities for the configuration x with and without the point u added. The Papangelou conditional intensity can be interpreted as follows: for any u ∈ R d and x ∈ N l f , λ (u, x)du corresponds to the conditional probability of observing a point in a ball of volume du around u given the rest of the point process is x.

Pairwise interaction point processes in R d can be defined and characterized through the Papangelou conditional intensity (see [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]) which is a function λ : R d ×N l f → R + . The Georgii-Nguyen-Zessin (GNZ) formula (see [START_REF] Papangelou | The Armenian Connection: Reminiscences from a time of interactions[END_REF], [START_REF] Zessin | Der Papangelou prozess[END_REF], [START_REF] Georgii | Canonical and grand canonical Gibbs states for continuum systems[END_REF], [START_REF] Nguyen | Integral and differential characterizations Gibbs processes[END_REF]) states that for any nonnegative measurable function

h on R d × N l f E ∑ u∈X h(u, X \ u) = E R d h(u, X)λ (u, X)du. (2.1)
Using induction we obtain the iterated GNZ-formula: for nonnegative func-

tions h : (R d ) n × N l f -→ R E = ∑ u 1 ,...,u n ∈X h(u 1 , . . . , u n , X \ {u 1 , . . . , u n }) = . . . E h(u 1 , . . . , u n , X)λ (u 1 , . . . , u n , X)du 1 . . . du n (2.2)
where λ (u 1 , . . . , u n , x) is Papangelou conditional intensity and is defined (not uniquely) by λ (u 1 , . . . , u n , x) = λ (u 1 , x)λ (u 2 , x ∪ {u 1 }) . . . λ (u n , x ∪ {u 1 , . . . , u n-1 }). Examples of Gibbs point process models and their conditional intensities are presented in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF], [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] and [START_REF] Møller | Modern statistics for spatial point processes[END_REF]. The formula (2.2) will also be used extensively later in this document.

Assumptions and the main results

For the general pairwise interaction process the conditional intensity is

λ (u, x) = g 0 (u) exp -∑ v∈x\u g 0 (u, v)
and note that g 0 (u, v) = g 0 (v, u) (i.e. symmetric pairwise interaction).

Alternatively, if g 0 (u) is a constant and g 0 (u, v) = g(vu) is translation invariant, then a pairwise interaction point process is called stationary.

Throughout this paper, we define a stationary pairwise interaction point process via the Papangelou conditional intensity at a location u given by

λ β (u, x) = β exp -∑ v∈x\u g(v -u) (3.1)
where β is the true value of the Poisson intensity parameter, g represents a pairwise interaction potential. Specifically, g is nonnegative measurable function; the practical significance of this assumption is that pairwise interaction processes are a wide and useful class of models for spatial point processes with inhibition between points; for discussion see [START_REF] Diggle | A nonparametric estimator for pairwise-interaction point processes[END_REF] and [START_REF] Stoyan | Stochastic geometry and its applications[END_REF]. We denote G = exp(-g) the pairwise interaction function. The basic assumption throughout this paper is the Papangelou conditional intensity (3.1) has a finite range R, i.e.

λ β (u, x) = λ β (u, x ∩ B(u, R)), (3.2)
for any u ∈ R d , x ∈ N l f , where B(u, R) is the closed ball in R d with center u and radius R.

Suppose that a single realization x of a point process X is observed in a bounded window Λ n ∈ B d 0 where (Λ n ) n≥1 is a sequence of cubes growing up to R d . We face a missing data problem, which in the spatial point process literature is referred to as a problem of edge effects, we can avoid this problem by reducing the window by introducing the 2R-interior of the cubes Λ n , i.e. Λ n,R = {u ∈ Λ n : B(u, 2R) ⊂ Λ n } and assume this has non-zero area. The choice of Λ n,R is induced by the fact that when u is on the edge of Λ n,R and v ∈ B(u, R), then we can observe the point v. Various edge corrections have been suggested by [START_REF] Ripley | Edge effects in spatial stochastic processes[END_REF], [START_REF] Biber | Ein verfahren zum ausgleish von rabdeffekten bei der berechnung von konkurrenzindizes[END_REF], [START_REF] Kelly | Statistical Inference for Spatial Processes[END_REF], [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] and [START_REF] Tukey | Spectral Analysis Time Series[END_REF]). We assume that the support of the interaction function G = exp(-g) is

T = {t ∈ R d ; g(t) > 0, for ||t|| < R}. Throughout in this paper, h is a nonnegative measurable function defined for all w ∈ R d , x ∈ N l f , by h(w, x) = 11(x ∩ B(w, R) = / 0), (3.3)
and we also introduce the following function

F(o, w) = E[h(o, X)h(w, X)] = P(X ∩ B(o, R) = / 0, X ∩ B(w, R) = / 0). (3.4)
To estimate the function β F(o,t), we propose empiric estimator Fn (t) defined for t ∈ T by

Fn (t) = 1 |Λ n,R | ∑ u∈X Λ n,R h(u, X \ {u})h(t + u, X \ {u}). (3.5)
To estimate the function β 2 G(t) F(o,t), we propose kernel-type estimator H n (t) defined for t ∈ T by

H n (t) = 1 b d n |Λ n,R | = ∑ u,v∈X v-u∈B(o,R) 11 Λ n,R (u)h(u, X \ {u, v})h(v, X \ {u, v})K v -u -t b n , (3.6
) where = over the summation sign means that the sum runs over all pairwise different points u; v in X and K : R d → R denotes a smoothing kernel function associated with a sequence (b n ) n≥1 of bandwidths. Plugging in the above estimators (3.5) and (3.6) and with the convention c/0 = 0 for all real c, we suggest a new edge-corrected nonparametric estimator G n (t) for β G(t) for t ∈ T by The kernel function K : R d -→ R is nonnegative and bounded with bounded support and satisfies:

G n (t) = H n (t) Fn (t) . ( 3
R d K(z)dz = 1, R d K 2 (z)dz < ∞. Condition K(d, m). Let z = (z 1 , . . . , z d ) , z i ∈ R, R d z i 1 1 . . . z i d d K(z 1 , . . . , z d )dz 1 . . . dz d = 0, for 0 < d ∑ j=1 i j < m,
and

R d |z| m K(z)dz < ∞.

Asymptotic behaviour of the variance and the covariance of the kernel-type estimator

In this section, we investigate the asymptotic behaviour of the variance and the covariance of the kernel-type estimator (3.6) estimating β 2 G(t) F(o,t). The latter is based on the knowledge of the iterated GNZ-formula (2.2).

Theorem 3.1. Consider a stationary pairwise interaction point process X in R d with Papangelou conditional intensity (3.1) satisfying condition (3.2). Further let the kernel function K satisfy Condition K(d). We have

lim n→∞ b d n |Λ n,R | Var H n (t) = β 2 G(t) F(o,t) R d K 2 (z)dz at any continuity point t ∈ T \{o} of G F.
The following theorem presents an asymptotic representation of the covariance of the kernel-type estimator (3.6) in two points t 1 ,t 2 in T . Theorem 3.2. Consider a stationary pairwise interaction point process X in R d with Papangelou conditional intensity (3.1) satisfying condition (3.2). Further let the kernel function K satisfy Condition K(d) and let t 1 = t 2 ∈ T . We have

lim n→∞ b d n |Λ n,R | Cov H n (t 1 ), H n (t 2 ) = 0.
The proof of Theorem 3.2 is similar to that of Theorem 3.1, and is omitted.

Asymptotic normality of a nonparametric estimator

Asymptotic normality of a nonparametric estimate of the pairwise interaction function G for a stationary pairwise interaction point process is essentially based on a multivariate central limit theorem for the estimator H n (t) of β 2 G(t) F(o,t). The proving idea of the below Theorem 3.3 consists in approximating the sequence H n (t) by a triangular array of m n -dependent random fields. The corresponding Lindeberg-type CLT for m n -dependent random fields has been proved in [START_REF] Wang | On the asymptotic normality of kernel density estimators for linear random fields[END_REF]. Now, we assume that the domain Λ n is divided into a fixed number of subdomains as follows Λ n = ∪ i∈I n Λ i , following [START_REF] Preston | Random fields[END_REF], [START_REF] Klein | Dobrushin uniqueness techniques and the decay of correlation in continuum statistical mechanics[END_REF], [START_REF] Jensen | On asymptotic normality of pseudo likelihood estimates of pairwise interaction processes[END_REF] we will describe a point process in R d as lattice process by means of this decomposition

Λ i = {ξ ∈ R d ; q(i j -1 2 ) ≤ ξ j ≤ q(i j + 1
2 ), j = 1, . . . , d} for a fixed number q > 0, i = (i 1 , . . . , i d ), and setting X i = X Λ i , i ∈ Z d , this becomes a Gibbs lattice field. We will consider estimation of

β 2 G(t) F(o,t) from H n (t), where the process is ob- served in Λ n,R = ∪ i∈ I n Λ i , where I n = {i ∈ I n ; |i -j| ≤ 1, for all j ∈ I n },
H n (t) = 1 b d n |Λ n,R | ∑ i∈ I n = ∑ u∈X i ,v∈X v-u∈B(o,R) h(u, X \ {u, v})h(v, X \ {u, v})K v -u -t b n .
We consider in this work the field {X i } i∈Z d (d ∈ N) in form of

X i = F(ε i-k , k ∈ Z d ), i ∈ Z d . (3.8) For each m ∈ N, i ∈ Z d , we write X i,m = F(ε i-k , k ∈ 0, m -1 d ),
X i = F(ε i-k , k ∈ Z d ),
where

ε j = ε j if j = 0 ε 0 if j = 0.
Let i ∈ Z d and p > 0 be fixed. We define δ i,p = X i -X i p , where • p is the usual L p -norm. Thus, [START_REF] Machkouri | Kernel density estimation for stationary random fields[END_REF] obtained this sufficient condition:

∑ i∈Z d |i| 5d/2 δ i,p < ∞ (3.9)
for a random field of the form (3.8) for the kernel density estimator to be asymptotically normal.

The desired sequence (m n ) n≥1 can be chosen as

m n = max v n , 1 b 3 n ∑ |i|>v n |i| 5d 2 δ i,2 1 3d
+ 1

where

v n = [b -1 2d
n ] and [•] denotes the integer part function. In order to establish the asymptotic normality of H n , we need additional assumptions:

Condition W . There exists a sequence of integers (m n ) n≥1 such that m n → ∞ as n → ∞, and the following limits hold Then, for any positive integer s and any distinct points t 1 , . . . ,t s in T \{o}, we have

m d n b d n → 0 and 1 (m n b n ) 3d/2 ∑ |i|>m n |i| 5d 2 δ i,2 → 0 if (3.
b d n |Λ n,R | H n (t i ) -E H n (t i ) s i=1 d -→ N (0, σ H )
where the covariance matrix σ H = [r i, j ] i, j=1,...,s is given by

r i,i = β 2 G(t i ) F(o,t i ) R d
K 2 (z)dz, for i = 1, . . . , s and r i, j = 0 if i = j.

d n |Λ n,R | H n (t i ) -β 2 G(t i ) F(o,t i ) s i=1 d -→ N (0, σ H )
where the covariance matrix σ H = [r i, j ] i, j=1,...,s is given by

r i,i = β 2 G(t i ) F(o,t i ) R d
K 2 (z)dz, for i = 1, . . . , s and r i, j = 0 if i = j.

The estimator (3.5) turns out to be unbiased estimator of β F(o,t) and strongly consistent (the uniform strong consistency) as n tends infinity, since a classical ergodic theorem for spatial point processes obtained in [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF]. This implies the following:

d n |Λ n,R | G n (t i ) -β G(t i ) s i=1 d -→ N (0, σ G )
where the covariance matrix σ G = [σ i, j ] i, j=1,...,s is given by Proof. We determine the asymptotic behaviour of the variance of the estimator H n (t). For this we use the following result.

σ i,i = G(t i ) F(o,t i ) R d K 2 (z)dz for i = 1, . . . , s and σ i, j = 0 if i = j.
Var = ∑ u,v∈X f (u, v, X\{u, v}) = E R 2d f (u, v, X) f (u, v, X) + f (v, u, X) λ (u, v, X)dudv + E R 3d f (u, v, X) f (v, w, X) + f (w, v, X)]λ (u, v, w, X)dudvdw + E R 3d f (u, v, X) f (u, w, X) + f (w, u, X) λ (u, v, w, X)dudvdw + E R 4d f (u, v, X) f (w, y, X)λ (u, v, w, y, X)dudvdwdy - R 4d E f (u, v, X)λ (u, v, X) E f (w, y, X)λ (w, y, X) dudvdwdy.
Proof. After a simple calculation, but prolix (see [START_REF] Jolivet | Upper bounds of the speed of convergence of moment density estimators for stationary point processes[END_REF] and [START_REF] Heinrich | Asymptotic Gaussianity of some estimators for reduced factorial moment measures and product densities of stationary Poisson cluster processes[END_REF]), we obtain the following decomposition:

= ∑ u,v∈X f (u, v, X\{u, v}) 2 = = ∑ u,v∈X f 2 (u, v, X\{u, v}) + = ∑ u,v∈X f (u, v, X\{u, v}) f (v, u, X\{u, v}) + = ∑ u,v,w∈X f (u, v, X\{u, v, w}) f (v, w, X\{u, v, w}) + = ∑ u,v,w∈X f (u, v, X\{u, v, w}) f (w, v, X\{u, v, w}) + = ∑ u,v,w∈X f (u, v, X\{u, v, w}) f (u, w, X\{u, v, w}) + = ∑ u,v,w∈X f (u, v, X\{u, v, w}) f (w, u, X\{u, v, w}) + = ∑ u,v,w,y∈X f (u, v, X\{u, v, w, y}) f (w, y, X\{u, v, w, y}). (4.1)
On the other hand using standard definition of variance, we have

Var = ∑ u,v∈X f (u, v, X\{u, v}) = E = ∑ u,v∈X f (u, v, X\{u, v}) 2 -E = ∑ u,v∈X f (u, v, X\{u, v}) 2 .
(4.2) After rearrangement with the formula of GNZ (2.2) and formulas (4.1) and (4.2), we obtain the desired result.

In the sequel, we denote

H(u 1 , u 2 , . . . , u s , X) = h(u 1 , X)h(u 2 , X) . . . h(u s , X), (4.3)
where h is given by (3.3),

I R (u 1 , u 2 , . . . , u s ) = 11 u 1 ∈ B(o, R), u 2 ∈ B(o, R), . . . , u s ∈ B(o, R) (4.4)
and we keep in mind the following function

F(u 1 , u 2 , . . . , u s ) = E[H(u 1 , u 2 , . . . , u s , X)] = E[h(u 1 , X)h(u 2 , X) . . . h(u s , X)]. (4.5) Substitute f (u, v, x) = 11 Λ n,R (u)I R (v -u)H(u, v, x)K v -u -t b n
the variance term of H n (t) is now expanded using Lemma 1, therefore we find that

Var H n (t) = T 1 + T 2 + T 3 + T 4 + T 5 + T 6 + T 7 + T 8 ,
with

T 1 = 1 b 2d n |Λ n,R | 2 E R 2d 11 Λ n,R (u)I R (v -u)H(u, v, X)K 2 v -u -t b n λ β (u, v, X)dudv, T 2 = 1 b 2d n |Λ n,R | 2 E Λ 2 n,R I R (v -u, u -v)H(u, v, X)K v -u -t b n K u -v -t b n λ β (u, v, X)dudv, T 3 = 1 b 2d n |Λ n,R | 2 E R d Λ 2 n,R I R (v -u, w -v)H(u, v, w, X)K v -u -t b n K w -v -t b n × λ β (u, v, w, X)dudvdw, T 4 = 1 b 2d n |Λ n,R | 2 E R d Λ 2 n,R I R (v -u, v -w)H(u, v, w, X)K v -u -t b n K v -w -t b n × λ β (u, v, w, X)dudvdw, T 5 = 1 b 2d n |Λ n,R | 2 E R 2d Λ n,R I R (v -u, w -u)H(u, v, w, X)K v -u -t b n K w -u -t b n × λ β (u, v, w, X)dudvdw, T 6 = 1 b 2d n |Λ n,R | 2 E R d Λ 2 n,R I R (v -u, u -w)H(u, v, w, X)K v -u -t b n K u -w -t b n × λ β (u, v, w, X)dudvdw, T 7 = 1 b 2d n |Λ n,R | 2 E R 4d 11 Λ n,R (u)11 Λ n,R (w)I R (v -u, y -w)H(u, v, w, y, X)K v -u -t b n × K y -w -t b n λ β (u, v, w, y, X)dudvdwdy
and

T 8 = - 1 b 2d n |Λ n,R | 2 R 4d 11 Λ n,R (u)11 Λ n,R (w)I R (v -u, y -w)K v -u -t b n K y -w -t b n × E H(u, v, X)λ β (u, v, X) E H(w, y, X)λ β (w, y, X) dudvdwdy.
The asymptotic behaviour of the leading term T 1 is obtained by applying the second order Papangelou conditional intensity given by:

λ β (u, v, x) = λ β (u, x)λ β (v, x ∪ {u}) for any u, v ∈ R d and x ∈ N l f .
And using the finite range property (3.2) for each function λ β (u, x) and λ β (v, x ∪ {u}). We recall when x = / 0, this implies that

λ β (u, / 0) = β and λ β (v, / 0 ∪ {u}) = β G(v -u) for all u, v ∈ R d .
By stationarity of X, it results

T 1 = 1 b 2d n |Λ n,R | 2 E R 2d 11 Λ n,R (u)I R (v -u)H(u, v, X)K 2 v -u -t b n λ β (u, v, X)dudv = β 2 b 2d n |Λ n,R | 2 R 2d 11 Λ n,R (u)I R (v -u) F(u, v)K 2 v -u -t b n G(v -u)dudv = β 2 b 2d n |Λ n,R | R d I R (s) F(o, s)K 2 s -t b n G(s)ds = β 2 b d n |Λ n,R | R d I R (b n z + t) F(o, b n z + t)K 2 (z)G(b n z + t)dz.
The continuity of G F in t and the boundedness conditions on the kernel function yield the desired result by dominated convergence theorem, i.e.

lim n→∞ b d n |Λ n,R |T 1 = β 2 G(t) F(o,t) R d K 2 (z)dz.
The following lemma will justify the applicability of dominated convergence theorem in the proofs of the next asymptotic results.

lim n→∞ Λ n ∩ (Λ n -b n y -z) |Λ n | = 1.
We will now show that all the other terms of the variance of the estimator H n (t) vanish. With similar arguments to those used to obtain the results above, we treat the term T 2 .

T 2 = 1 b 2d n |Λ n,R | 2 E R 2d 11 Λ n,R (u)11 Λ n,R (v)I R (v -u, u -v)H(u, v, X)K v -u -t b n × K u -v -t b n λ β (u, v, X)dudv = β 2 b 2d n |Λ n,R | 2 R 2d 11 Λ n,R (u)11 Λ n,R (v)I R (v -u, u -v) F(o, v -u)K v -u -t b n × K u -v -t b n G(v -u)dudv = β 2 b 2d n |Λ n,R | R d |Λ n,R ∩ (Λ n,R -s)| |Λ n,R | I R (s, -s) F(o, s)K s -t b n K -s -t b n G(s)ds = β 2 b d n |Λ n,R | R d |Λ n,R ∩ (Λ n,R -(b n z + t))| |Λ n,R | I R (b n z + t, -b n z -t) F(o, b n z + t)K(z) × K(-z - 2t b n )G(b n z + t)dz.
The boundedness conditions on the kernel and under condition t = o imply that b d n |Λ n,R |T 2 -→ 0 as n → ∞, by dominated convergence theorem and the geometric properties of Λ n,R (Lemma 2).

For the asymptotic behaviour of the third term T 3 , we remember the third order Papangelou conditional intensity by

λ β (u, v, w, x) = λ β (u, x)λ β (v, x ∪ {u})λ β (w, x ∪ {u, v})
for any u, v, w ∈ R d and x ∈ N l f . Using the finite range property (3.2) for each function λ β (u, x), λ β (v, x ∪ {u}) and λ β (w, x ∪ {u, v}), i.e.

λ β (u, X) = λ β (u, X ∩ B(u, R)) = β when X ∩ B(u, R) = / 0, λ β (v, X ∪ {u}) = λ β (v, (X ∪ {u}) ∩ B(v, R)) = β G(v -u) when X ∩ B(v, R) = / 0, v -u ∈ B(o, R),
and when X ∩ B(w, R) = / 0 and wv ∈ B(o, R), we have

λ β (w, X ∪ {u, v}) = λ β (w, (X ∪ {u, v}) ∩ B(w, R)) = λ β (w, v, {u} ∩ B(w, R)).
Since X is a point process to interact in pairs, the interaction terms due to triplets or higher order are equal to one, i.e. G(y) = 1 when n(y) > 2, for / 0 = y ⊆ x. Then,

λ β (w, v, {u} ∩ B(w, R)) = β G(w -v) if u / ∈ B(w, R), β otherwise and λ β (u, v, w, / 0) = β 3 G(v -u)G(w -v) if u / ∈ B(w, R) β 3 G(v -u) otherwise. (4.6)
Which ensures that λ β (u, v, w, / 0) is a function that depends only variables vu, wv, denoted by G 3 (vu, wv). In this way from the definition of H (resp. I R ) given by (4.3) (resp.(4.4)) and by stationarity of X, it follows that

T 3 = 1 b 2d n |Λ n,R | R 2d |Λ n,R ∩ (Λ n,R -v)| |Λ n,R | Ḡ3 (v, w -v)K v -t b n K w -v -t b n dvdw = 1 |Λ n,R | R 2d |Λ n,R ∩ (Λ n,R -(b n z + t))| |Λ n,R | Ḡ3 (b n z + t, b n z + t)K(z)K(z )dzdz , where Ḡ3 (v, w -v) = I R (v, w -v) F(-v, w -v)G 3 (v, w -v)
, where F is given by (4.5). By dominated convergence theorem and using the Lemma 2, we get b d n |Λ n,R |T 3 -→ 0 as n → ∞. Analogously, one may show the other terms T 4 , T 5 and T 6 .

For the asymptotic behaviour of the leading term T 7 , let u, v, w, y any points in R d and x ∈ N l f , we define Papangelou conditional intensity of the fourth order by λ β (u, v, w, y, x) = λ β (u, x)λ β (v, x ∪ {u})λ β (w, x ∪ {u, v})λ β (y, x ∪ {u, v, w}).

Next we introduce the finite range property (3.2) and reasoning analogous with the foregoing on λ β (u, v, w, / 0) is defined by (4.6), which ensures that λ β (u, v, w, y, / 0) is a function that depends only variables vu, yw, wu, wv, denoted by G 4 (vu, yw, wu, wv). Further by the stationarity of X, this implies that

T 7 = 1 b 2d n |Λ n,R | 2 E R 4d 11 Λ n,R (u)11 Λ n,R (w)I R (v -u, y -w)H(u, v, w, y, X)K v -u -t b n × K y -w -t b n λ β (u, v, w, y, X)dudvdwdy = 1 b 2d n |Λ n,R | R 3d |Λ n,R ∩ (Λ n,R -w)| |Λ n,R | Ḡ4 (v, y -w, w, w -v)K v -t b n K y -w -t b n dvdwdy = 1 |Λ n,R | R 3d |Λ n,R ∩ (Λ n,R -w)| |Λ n,R | Ḡ4 (b n z + t, b n z + t, w, w -b n z -t)K(z)K(z )dzdz dw, where Ḡ4 (v, y -w, w, w -v) = I R (v, y -w) F(v, w, y)G 4 (v, y -w, w, w -v). Similar arguments show that b d n |Λ n,R |T 7 -→ 0 as n → ∞.
We now finish the proof of the asymptotic behaviour of the term T 8 . Under the conditions imposed in Theorem 3.1, we have

T 8 = - 1 b 2d n |Λ n,R | 2 R 4d 11 Λ n,R (u)11 Λ n,R (w)I R (v -u, y -w)K v -u -t b n × K y -w -t b n E H(u, v, X)λ β (u, v, X) E H(w, y, X)λ β (w, y, X) dudvdwdy = - 1 b 2d n |Λ n,R | 2 R 4d 11 Λ n,R (u)11 Λ n,R (w)K v -u -t b n K y -w -t b n × Ḡ5 (u, v, w, y, v -u, y -w)dudvdwdy = - 1 b 2d n |Λ n,R | × R 3d |Λ n,R ∩ (Λ n,R -w)| |Λ n,R | Ḡ5 (o, v, w, y, v, y -w)K v -t b n K y -w -t b n dvdwdy = - 1 |Λ n,R | R 3d |Λ n,R ∩ (Λ n,R -w)| |Λ n,R | Ḡ5 (o, b n z + t, w, b n z + t + w, b n z + t, b n z + t) × K(z)K(z )dwdzdz , where Ḡ5 (u, v, w, y, v -u, y -w) = I R (v -u, y -w) F(u, v) F(w, y)G(v -u)G(y -w)
, where F is given by (4.5). By dominated convergence theorem, we get the desired result.

Proof of Theorem 3.3

Proof. To prove the Theorem 3.3, we need (using the well-known Cramér-Wold (1936)). Without loss of generality, we consider only the case s = 2. Let λ 1 and λ 2 be two constants with λ 1 + λ 2 = 0 and let the points t 1 ,t 2 ∈ T \{o} be fixed, t 1 = t 2 and define the triangular array of random variables 1

|Λ n,R | 1/2 ∑ i∈ I n Ȳn,i = λ 1 (b d n |Λ n,R |) 1/2 H n (t 1 ) -E[ H n (t 1 )] + λ 2 (b d n |Λ n,R |) 1/2 H n (t 2 ) -E[ H n (t 2 )]
where

Ȳn,i = λ 1 Y i (t 1 ) -E[Y i (t 1 )] + λ 2 Y i (t 2 ) -E[Y i (t 2 )]
and for all j = 1, 2, we consider

Y i (t j ) = (b d n ) -1/2 = ∑ u∈X i ,v∈X v-u∈B(o,R) h(u, X \ {u, v})h(v, X \ {u, v})K v -u -t j b n .
The proving idea the below Theorem 3.3 consists in approximating the sequence Ȳn,i by a triangular array of m-dependent random fields. The corresponding Lindebergtype CLT m-dependent random fields has been proved by [START_REF] Wang | On the asymptotic normality of kernel density estimators for linear random fields[END_REF]. Therefore, it proves convenient to switch notation from the text and to define

Y n,i = λ 1 Y i (t 1 ) -E[ Y i (t 1 )] + λ 2 Y i (t 2 ) -E[ Y i (t 2 )]
and for all j = 1, 2, we consider

Y i (t j ) = (b d n ) -1/2 = ∑ u∈X i,

Proof of Proposition 1

Proof. Taking into account of { Y n,i } i∈Z d are m n -dependent random fields, we apply limit theorem for stationary triangular arrays of m-dependent random fields (Theorem 2 established by [START_REF] Wang | On the asymptotic normality of kernel density estimators for linear random fields[END_REF]). We need to verify that the two fundamental conditions (4.9) and (4.10), i.e. it suffices to show, for some sequence of increasing integers {l n } n∈N such that m n /l n → 0 and l n /n → 0 as n → ∞,

Let B d

  be the Borel σ -algebra (generated by open sets) in R d (the d-dimensional space) and B d

. 7 )

 7 Moreover, we have to impose certain natural restrictions on the kernel function K and the sequence (b n ) n≥1 . Condition K(d). The sequence of bandwidths (b n ) n≥1 , is a sequence of positive real numbers satisfying: lim n→∞ b n = 0 and lim n→∞ b d n |Λ n,R | = ∞.

  and the norm is | j| = max{| j 1 |, . . . , | j d |} and assume that I n increases towards Z d and write

Theorem 3 . 3 .

 33 Consider a stationary pairwise interaction point process X in R d with Papangelou conditional intensity (3.1) satisfying condition (3.2). Further let the kernel function K satisfy Condition K(d) and assume that Condition W holds.

Corollary 1 .

 1 Let, in addition to the assumptions of Theorem 3.3 and the kernel function K satisfy Condition K(d, m) and let b d+2m n |Λ n,R | → 0 as n → +∞. Furthermore if G(t) F(o,t) has bounded and continuous partial derivatives of order m in some neighborhood of the points t 1 , . . . ,t s . Then b

Theorem 3 . 4 .

 34 Consider a stationary pairwise interaction point process X in R d with Papangelou conditional intensity (3.1) satisfying condition (3.2). Further let the kernel function K satisfy Condition K(d), Condition K(d, m) and assume that Condition W and b d+2m n |Λ n,R | → 0 as n → ∞ hold. Then for any positive integer s and any distinct points t 1 , . . . ,t s in T \{o}, b

Lemma 1 .

 1 Consider any Gibbs point process X in R d with Papangelou conditional intensity λ . Let f : R d × R d × N l f -→ R be a nonnegative and measurable function. Then

Lemma 2 .

 2 ([8], Lemma A.2) Let (Λ n ) n∈N be a sequence of convex sets in R d growing up to R d and let y, z ∈ R d . Let (b n ) n∈N be a sequence of positive real numbers with b n -→ n→∞ 0. Then we have

  i > εn d/2 = 0 for all ε > 0 (4.10) where a, b ≡ {a, a + 1, . . . , b} for a, b ∈ Z. Taking into account the definition of ∑ i∈ 1,l n d Y n,i we get lim

  where {ε i } i∈Z d are i.i.d. random variables and F is measurable function. Let {ε i } i∈Z d be an i.i.d. copy of {ε i } i∈Z d and consider for all positive integer n the coupled version X i defined by

  m n ,v∈X v-u∈B(o,R) h(u, X \ {u, v})h(v, X \ {u, v})K vut j b n . By construction, { Y n,i } i∈Z d are m n -dependent, in the sense that Y n,i and Y n, j are independent if |i -j| ∞ ≥ m n , where | j| ∞ = max{| j 1 |, . . . , | j d |}. We will use { Y n,i ; i ∈ Z d } n∈N to approximate { Ȳn,i ; i ∈ Z d } n∈N . We decompose 1 |Λ n,R | 1/2 ∑ Ȳn,i -Y n,i .Since K and h are bounded and applying Proposition 1 obtained[START_REF] Machkouri | A central limit theorem for stationary random fields[END_REF] and under (3.10) of Condition W we know that 1|Λ n | 1/2 ∑Then, it suffices to establish the following result. Proposition 1. Under (3.11) of Condition W . As n → ∞, we have 1 |Λ n,R | 1/2∑ = β 2 λ 2 1 F(o,t 1 )G(t 1 ) + β 2 λ 2 2 F(o,t 2 )G(t 2 )

	i∈ I n	Ȳn,i -Y n,i	2	≤ c	|λ 1 | + |λ 2 | (m n b n ) 3d/2 ∑ |i|>m n	|i|	5d 2 δ i,2 = o(1).
					Y n,i	d -→ N 0, σ 2 ,
				i∈ I n
	where
	σ 2 R d	K 2 (z)dz.	(4.7)
	i∈ I n	Ȳn,i =	1 |Λ n,R | 1/2 ∑ i∈ I n	Y n,i +	1 |Λ n,R | 1/2 ∑ i∈ I n

  d n l d n λ 2 1 Var( H l n (t 1 )) + λ 2 2 Var( H l n (t 2 )) + 2λ 1 λ 2 Cov( H l n (t 1 ), H l n (t 2 )) .

Acknowledgements: The research was supported by laboratory Jean Kuntzmann, Grenoble University, France.

where σ 2 is defined by (4.7).

To prove (4.10), we write for all n ∈ N. For p > 2, observe that, since K and h are bounded,

.

By the above inequalities, we have obtained

. Now, (3.11) of Condition W entails that l n can be chosen so that in addition to (4.8),

This completes the proof of Theorem 3.3.

Proof of Corollary 1

Proof. Let t ∈ {t 1 , . . . ,t s } be fixed. Using the second-order Georgii-Nguyen-Zessin formula (2.2), and using the finite range property (3.2) and by stationarity of X, we obtain

By Taylor expansion of the integrand in neighborhood of t and making use of Condition K(d, m), Condition K(d) and the function G(t) F(o,t) has bounded and continuous partial derivatives of order m of in some neighborhood of the point t, we get the following rate of convergence

By the application of Theorem 3.3 we conclude the result announced.

Proof of Theorem 3.4

Proof. Let t ∈ {t 1 , . . . ,t s } be fixed. We may split the difference G n (t)β G(t) as follows:

n .

Hence, from [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF] ergodic theorem, and additionally assume that P is ergodic, Fn (t) converges almost surely to β F(o,t). Note that there exists at least one stationary Gibbs measure. If this measure is unique, it is ergodic. Otherwise, it can be represented as a mixture of ergodic measures (see [START_REF] Georgii | Gibbs measures and Phase Transitions[END_REF], Theorem 14.10). Therefore, we can assume for this proof, that P is ergodic. By Theorem 3.3, we have that as

n tends to a Gaussian distribution with zero mean. From strong consistency of the estimator Fn (t), and by inserting (4. which completes the proof.


Nadia Morsli
email: nadia.morsli@imag.fr

Uniform strong consistency of nonparametric estimator for stationary pairwise interaction point processes under mixing conditions

Keywords: kernel-type estimator, mixing random fields, pairwise interaction spatial point process, Papangelou conditional intensity, rates of converge, uniform strong consistency

We suggest a new nonparametric estimate of the interaction function of pairwise interaction point process on R d such that its Papangelou conditional intensity is translation invariant and satisfied finite-range interaction. We prove uniform strong consistency of the nonparametric estimate of the interaction function . Sufficient conditions of the rates of uniform strong consistency for the resulting kernel-type estimator are valid by the nonuniform mixing condition for a stationary point process, which can be viewed interchangeably as a lattice field in Z d .

Introduction

Point processes are well studied objects in probability theory and statistics, which is of interest in such diverse disciplines as forestry, seismology, spatial epidemiology, materials science, astronomy, geography, ecology and other disciplines. In this work we will develop some recent results of nonparametric statistical spatial processes. Functional data analysis is a typical issue in modern statistics. During the last years, many papers have been devoted to theoretical results or applied studies on models involving functional data. The main goal of this paper is to present a nonparametric estimation of interaction function for stationary pairwise interaction point processes (they have been introduced in statistical literature by Ripley and Kelly [START_REF] Ripley | Markov point processes[END_REF], Daley and Vere-Jones [START_REF] Daley | An introduction to the Theory of Point Processes[END_REF] and Georgii [START_REF] Georgii | Gibbs measures and Phase Transitions[END_REF]) defined in the space R d for which its Papangelou conditional intensity (Papangelou [17]). In general it is not possible to deal with densities of infinite point processes, instead, the Papangelou conditional intensity becomes the appropriate starting point for modelling. A formal definition and interpretation of the Papangelou conditional intensity are given below. The estimate of the interaction function plays a very important role in the statistical analysis of spatial point pattern data. We explore the spatial dependence of random fields (mixing random fields) and we discuss the uniform strong consistent estimation. Many modelizations of the mixing random fields have been proposed through in literature, we refer to Bradley [START_REF] Bradley | Introduction to Strong Mixing Conditions[END_REF], Rosenblatt [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF], I.A.Ibragimov [START_REF] Ibragimov | Some limit theorems for stationary processes[END_REF]. Bosq [START_REF] Bosq | La statistique non paramétrique des processus[END_REF], Delecroix [START_REF] Delecroix | [END_REF], Heinrich and Liebscher [START_REF] Heinrich | Strong convergence of kernel estimators for product densities of absolutely regular point processes[END_REF], and Pelidrad [START_REF] Pelidrad | Properties of uniform consistency of the kernel estimators of density and of regression functions under dependence assumptions[END_REF]. The mixing hypothesis was introduced to extend the central limit theorem for more general classes dependent random variables (see Ibragimov and Linnik [START_REF] Ibragimov | Independent and Stationary Sequences of Random Variables[END_REF]). It also helps to generalize almost sure convergence obtained with independent data.

The paper is organized as follows. In Section 2, we set up the generic notation. In Section 3, we present the model that has been considered for stationary pairwise interaction point process and we suggest a new nonparametric estimator of interaction function. We present some conditions that will be helpful in the following paper and we present the main results, We suggest a new nonparametric estimate of the interaction function of pairwise interaction point process on R d in the Papangelou conditional intensity and we prove uniform strong consistency of the resulting estimator. The Section 4 is devoted to the proofs. Our methods of proofs are based on new Kahane-Khintchine inequalities in Orlicz spaces induced by exponential Young functions for stationary uniformly strong mixing random fields, obtained by (El Machkouri [START_REF] El Machkouri | Kahane-Khintchine inequalities and functional central limit theorem for stationary real random fields[END_REF].

Generic notation and Basic tools

We introduce some necessary definitions and notations and recall some basic facts from the theory of Gibbs point processes. Let B d be the Borel σ -field on R d of general dimension d ≥ 2 and The Papangelou conditional intensity of Gibbs point processes X in R d (Møller and Waagepetersen [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]) is a function λ : R d ×N l f → R + and characterized by the Georgii-Nguyen-Zessin (GNZ) formula (see Papangelou [START_REF] Papangelou | The Armenian Connection: Reminiscences from a time of interactions[END_REF] and Zessin [START_REF] Zessin | Der Papangelou prozess[END_REF] for historical comments and Georgii [START_REF] Georgii | Canonical and grand canonical Gibbs states for continuum systems[END_REF] or Nguyen and Zessin [START_REF] Nguyen | Integral and differential characterizations Gibbs processes[END_REF] for a general presentation). The GNZ formula states that for any nonnegative measurable function

B d O = {Λ ∈ B d : Λ bounded} . For a set Λ ∈ B d , |Λ| is the positive volume of Λ. We define the space of locally finite point configurations in R d as N l f = {x ⊆ R d ; n(x ∩ Λ) < ∞, ∀Λ ∈ B d 0 } equipped with the smallest σ -field N l f containing all set of the form {x ∈ N l f : n(x ∩ Λ) = k} for k ∈ N 0 = {0, 1, 2,
h on R d × N l f E ∑ ξ ∈X h(ξ , X \ ξ ) = E R d h(ξ , X)λ (ξ , X)dξ . (2.1)
Heuristically, λ (ξ , x)dξ can be interpreted as the conditional probability of X having a point in an infinitesimal small region containing ξ and of size dξ given the rest of X is x. Using induction we obtain the second-order GNZ-formula

E = ∑ ξ ,η,∈X h(ξ , η, X \ {ξ , η}) = R d R d E h(ξ , η, X)λ (ξ , η, X)dξ dη (2.2) where λ (ξ , η, x) = λ (ξ , x)λ (η, x ∪ {ξ }), for ξ , η ∈ R d , x ∈ N l f with a nonnega- tive measurable function h(ξ , η, x) on R d × R d × N l f .
3 Main results

The model and its estimator

Throughout this paper, we shall focus on stationary pairwise interaction point processes (Daley and Vere-Jones [START_REF] Daley | An introduction to the Theory of Point Processes[END_REF], Georgii [START_REF] Georgii | Gibbs measures and Phase Transitions[END_REF]). We specify that the Papangelou conditional intensity is translation invariant for a configuration x ∈ N l f and a location ξ ∈ R d should be of the form

λ (ξ , x) = β exp -∑ η∈x\ξ V 2 (η -ξ ) (3.3)
where β is the true value of the intensity parameter. The function V 2 represents the pair potential defined on R d . To do this, we will assume the following condition on the potential

V 2 : 0 < V 2 (ξ ) < ∞. (3.4)
Condition 3.4 says that the interaction is purely repulsive. In physics applications, V 2 is interpreted as the energy of the configuration x. V = exp(-V 2 ) represents the pairwise interaction function. The Papangelou conditional intensity is said to have finite interaction range R with 0

< R < ∞, if λ (ξ , x) = λ (ξ , x ∩ B(ξ , R)), (3.5)
for any ξ ∈ R d , x ∈ N l f , where B(ξ , R) is the closed ball in R d with center ξ and radius R.

Our objective is to study problems of estimation of the interaction function V in the model (3.3). For this, let Λ n ∈ B d O , for n ∈ N, be a sequence, denotes a bounded observation window in such a way that

|Λ n | → ∞, as n → ∞. In most application examples Λ n cubes growing up to R d . Throughout in this paper, f is a nonnegative measurable function defined for all w, w ∈ R d , x ∈ N l f , by f (w, w , x) = 11(x ∩ B(w, R) = / 0, x ∩ B(w , R) = / 0)
and we also introduce the following function

C(o,t) = E[ f (o,t, X)] = P(X ∩ B(o, R) = / 0, X ∩ B(t, R) = / 0). (3.6)
Now, we suggest a new edge-corrected nonparametric estimator for the function β V (t) by

V n (t) = J n (t) C n (t) , for t ∈ R d (3.7)
where

J n (t) = = ∑ ξ ,η∈X∩Λ n 11(η -ξ ∈ B(o, R)) f (ξ , η, X \ {ξ , η})K η-ξ -t b n b d n |Λ n ∩ (Λ n + (η -ξ ))| (3.
C n (t) = ∑ ξ ∈X∩Λ n f (ξ ,t + ξ , X \ {ξ }) |Λ n ∩ (Λ n + ξ )| . (3.9)
In this case, the denominator equals 0, we choose (by convention) an usual way to define our estimator : V n (t) = 0.

Uniform strong consistency

The proof of uniform strong consistency of the nonparametric estimator V n (t) defined by (3.7), requires more delicate conditions. We have to impose certain restriction on the kernel function K and the sequence of bandwidths (b n ) n≥1 . The following assumptions are required.

(A1). The kernel function K : R d -→ [0, ∞) is bounded with bounded support, and satisfies:

Let ξ = (ξ 1 , . . . , ξ d ) , ξ i ∈ R, R d ξ i 1 1 . . . ξ i d d K( x i 1 , . . . , ξ d )dξ 1 . . . dξ d = 1 if i 1 = i 2 = . . . i d = 0, 0 if 0 < i 1 + i 2 + . . . + i d < s. R d ||ξ || i |K(ξ )|dξ < ∞ for i = 0 and i = s, R d K(ξ )dξ = 1.
(A2). The kernel function K satisfies a Lipschitz condition, i.e. there exists a constant L > 0 such that These assumptions are more than enough to guarantee rates of convergence of mean and a uniform strong consistency of the kernel-type estimator J n (t) (defined by (3.8)). Now, we establish thereafter interesting results for establishing the uniform strong consistency of the estimator V n (t).

K(ξ ) -K(η) ≤ L||ξ -η|| for any ξ , η ∈ R d . ( A3
Theorem 1. Let X be a stationary pairwise interaction point process in R d with Papangelou conditional intensity (3.3) satisfying conditions (3.4) and (3.5). Let the kernel function K satisfy (A1) and (A2). Let the bandwidths (b n ) n≥1 satisfy (A3). Furthermore, we assume that V (t) C(o,t) has bounded and continuous partial derivatives of order s on for any fixed compact set K 0 ⊂ R d . Then,

sup t∈K 0 V n (t) -β V (t) a.s -→ n→∞ 0.
The proof of uniform strong consistency of the nonparametric estimator V n (t) (Theorem 1) is based on some asymptotic properties of the estimators C n (t) and J n (t). For the estimator J n (t), we are exploring the spatial dependence of point processes, i.e. φ -mixing conditions.

Rates of uniform strong consistency of the kernel-type estimator J n (t)

In this section we discuss sufficient conditions which ensure uniform strong consistency (including a convergence rate) of the kernel-type estimator J n (t) under uniformly strong mixing (i.e. φ -mixing) random fields. We use Kahane-Khintchine inequalities in Orlicz spaces induced by exponential Young functions for stationary φ -mixing real random fields which satisfy some finite exponential moment condition. Now let us the nonuniform φ -mixing condition. Given two sub-σ -algebras U and V of F . The φ -mixing coefficients had been introduced by I.A.Ibragimov [START_REF] Ibragimov | Some limit theorems for stationary processes[END_REF] and can be defined by:

φ (U , V ) = sup{ P(V |U ) -P(V ) ∞ , V ∈ V }.
In the sequel, we consider the nonuniform φ 1,∞ (n)-mixing coefficients φ 1,∞ (n) defined for each positive integer n by

φ 1,∞ (n) = sup {φ (σ (X k ), F Γ ), k ∈ Z d , Γ ⊂ Z d , d(Γ, {k}) ≥ n}, where F Γ = σ (X i , i ∈ Γ) is the σ -algebra genereted by X i , i ∈ Γ, |Γ| is the cardinal of Γ and d(Γ 1 , Γ 2 ) = min{|i -k| : i ∈ Γ 1 , k ∈ Γ 2 }.
We say that the random field

(X i ) i∈Z d is φ -mixing if lim n→∞ φ 1,∞ (n) = 0.
As shown by (Dobrushin [6]; Georgii [START_REF] Georgii | Gibbs measures and Phase Transitions[END_REF]; Guyon [START_REF] Guyon | Random Fields on a Network: Modeling, Statistics and Applications[END_REF]), the weak dependence conditions based on the above mixing coefficients are satisfied by large classes of random fields including Gibbs fields. Now, we assume that the domain Λ n is divided into a fixed number of subdomains as follows Λ n = ∪ i∈Γ n Λ i , where Λ i be the unit cube centered at i ∈ Z d and assume that Γ n ⊂ Z d , such that |Γ n | → ∞, as n → ∞. We will consider estimation of β 2 V (t) C(o,t) from J n (t), where the process is observed in ∪ i∈Γ n Λ i as

J n (t) = ∑ i∈Γ n = ∑ ξ ,η∈X∩Λ i 11(η -ξ ∈ B(o, R)) f (ξ , η, X \ {ξ , η})K η-ξ -t b n b d n |Λ n ∩ (Λ n + (η -ξ ))| . Put J n,i = = ∑ ξ ,η∈X∩Λ i 11(η -ξ ∈ B(o, R)) f (ξ , η, X \ {ξ , η})K η -ξ -t b n
and we note that: Jn,i = J n,i -E(J n,i ).

Recall that a Young function ψ is a real convex nondecreasing function defined on R + which satisfies lim t→∞ ψ(t) = +∞ and ψ(0) = 0. We define the Orlicz space L ψ as the space of real random variables Z defined on the probability space (Ω, F , P) such that E[ψ(|Z|/c)] < +∞ for some c > 0. The Orlicz space L ψ equipped with the so-called Luxemburg norm . ψ defined for any real random variable Z by

Z ψ = inf{ c > 0 ; E[ψ(|Z|/c)] ≤ 1 }
is a Banach space. Let θ > 0. We denote by ψ θ the exponential Young function defined for any x ∈ R + by

ψ θ (x) = exp((x+ξ θ ) θ )-exp(ξ θ θ ) where ξ θ = ((1-θ )/θ ) 1/θ 11{0 < θ < 1}.
(3.10) On the lattice Z d we define the lexicographic order as follows: if i = (i 1 , . . . , i d ) and j = ( j 1 , . . . , j d ) are distinct elements of Z d , the notation i < lex j means that either i 1 < j 1 or for some p in {2, 3, . . . , d}, i p < j p and i q = j q for 1 ≤ q < p. Let the sets V 1 i ; i ∈ Z d be defined as follows:

V 1 i = { j ∈ Z d ; j < lex i}.
To establish the proof of Theorem 2, we need the following result.

Corollary 1. (El Machkouri [START_REF] El Machkouri | Kahane-Khintchine inequalities and functional central limit theorem for stationary real random fields[END_REF]). Let (X i ) i∈Z d be a zero mean stationary real random field which satisfies the assumption:

∃q ∈]0, 2[ ∃c > 0 E[exp(c|X 0 | θ (q) )] < +∞,
where θ (q) = 2q/(2q). There exists a positive universal constant M 1 (q) > 0 depending only on q such that for any family (a i ) i∈Z d of real numbers and any finite subset

Γ in Z d , ∑ i∈Γ a i X i ψ q ≤ M 1 (q) X 0 ψ θ (q) ∑ i∈Γ |a i | b i,q (X) 1/2 where b i,q (X) = |a i | +C(q) ∑ k∈V 1 0 |a k+i | φ ∞,1 (|k|).
C(q) is the positive constant depending only on q.

Theorem 2. Let the kernel function K satisfy (A2). Then, for any fixed compact set K 0 ⊂ R d , we have: If there exists 0 < q < 2 such that Jn,0 ∈ L ψ θ (q) and ∑

k∈Z d φ 1,∞ (|k|) < ∞. (3.11)
Then

sup t∈K 0 J n (t) -E( J n (t)) = O a.s. (log n) 1/q (b n √ n) d as n → ∞.
Corollary 2. If in addition the function V (t) C(o,t) has bounded and continuous partial derivatives of order s on K 0 ⊂ R d and we suppose that (A1) and (A3) are fulfilled. If there exists 0 < q < 2 such that Jn,0 ∈ L ψ θ (q) and ∑

k∈Z d φ 1,∞ (|k|) < ∞.
Then

sup t∈K 0 J n (t) -β 2 V (t) C(o,t) = O a.s. (log n) 1/q (b n √ n) d + O(b s n ) as n → ∞.

Proofs

Proof of Theorem 2. Let (a n ) n≥1 be sequence of positive numbers going to zero. Following Carbon et al. [START_REF] Carbon | Kernel density estimation for random fields[END_REF], the compact set K 0 can be covered by v n cubes T k having sides of length L n = a n b d+1 n and center at c k . Clearly there exists c > 0, such that v n ≤ c/L d n . We use the following classical decomposition sup

t∈K 0 | J n (t) -E J n (t)| ≤ max 1≤k≤v n sup t∈T k J n (t) -J n (c k ) + max 1≤k≤v n sup t∈T k E J n (t) -E J n (c k ) + max 1≤k≤v n J n (c k ) -E J n (c k ) = B 1 + B 2 + B 3 .
For any t ∈ T k , by assumption (A2), we derive that there exists constant L > 0 such that n sufficiently large

J n (t) -J n (c k ) ≤ La n 1 cn d = ∑ ξ ,η∈X∩Λ n 11(η -ξ ∈ B(o, R)) f (ξ , η, X \ {ξ , η}).
Using the spatial ergodic theorem of Nguyen and Zessin [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF], and by the secondorder Georgii-Nguyen-Zessin formula (2.2) with

h(ξ , η, X) = 11(η -ξ ∈ B(o, R)) f (ξ , η, X) it easily follows B 1 = O a.s. (a n ).
From L 1 -version of the ergodic theorem of Nguyen and Zessin [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF], it follows

B 2 = O a.s. (a n ).
Now, we focus on the term B 3 . In the sequel, the letter C denotes any generic positive constant. Let ε > 0 and t ∈ K 0 be fixed. We use the Tchebychev-Markov inequality, we have

P | J n (t) -E J n (t)| > εa n = P | ∑ i∈Γ n Jn,i | > εa n (nb n ) d ≤ exp - ε a n (nb n ) d || ∑ i∈Γ n Jn,i || ψ q + l q q E exp | ∑ i∈Γ n Jn,i | || ∑ i∈Γ n Jn,i || ψ q + l q q
. (4.12)

Therefore, we assume that there exists a real 0 < q < 2, such that Jn,0 ∈ L ψ θ (q) . As a direct application of Kahane-Khintchine inequality of Corollary 1 to the zero mean sequence Jn,i , for i ∈ Γ n , we obtain that:

P | ∑ i∈Γ n Jn,i | > εa n (b n n) d ≤ (1 + e l q q ) exp -   ε a n (b n n) d M 1 (q) Jn,0 2 ψ θ (q) (∑ i∈Γ n b i,q ( J)) 1/2 + l q   q where b i,q ( J) = 1 +C(q) ∑ k∈V 1 0 φ 1,∞ (|k|),
M 1 (q) and C(q) are positive universal constants depending only on q. Thus, under condition (3.11) and from the stationary of X and by definition of the norm || • || ψ θ (q) , we infer that there exists a constant C > 0 such that

P(| J n (t) -E J n (t)| > εa n ) ≤ (1 + e l q q ) exp - ε q a q n (b n √ n) dq C q .
Inserting a n = (log n) 1/q /(b n √ n) d , we see after a short calculation that

P(| J n (t) -E J n (t)| > εa n ) ≤ (1 + e l q q ) exp - ε q log n C q .
From the last equality, we find that sup

t∈K 0 P(| J n (t) -E J n (t)| > εa n ) ≤ (1 + e l q q ) exp - ε q log n C q . ( 4.13)
Since

P (|B 3 | > εa n ) ≤ v n sup t∈K 0 P(| J n (t) -E J n (t)| > εa n ),
using (4.13), it follows with Borel-Cantelli's lemma

P(lim sup n→∞ |B 3 | > εa n ) = 0
We conclude the proof of Theorem 2.

By Taylor expansion of the integrand in neighborhood of t and making use of (A1), (A3) and the function V (t) C(o,t) has bounded and continuous partial derivatives of order s of in B o (t, δ ) (for some δ > 0 ), we get the following rate of convergence

E J n (t) -β 2 V (t) C(o,t) = O(b s n ) as n → ∞.
By the Theorem 2, we complete the proof.

Proof of Theorem 1. We consider the decomposition

V n (t) -β V (t) = J n (t) -β 2 V (t) C(o,t) C n (t) + β V (t) C n (t) -β C(o,t) C n (t) . (4.14)
Now, we prove by the ergodic theorem (Nguyen and Zessin [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF]) that

C n (t) -β C(o,t)
a.s.

-→ 0, as n → ∞. (4.15)

Let (τ y ) y∈R d be the shift group, where τ y is the translation by the vector -y ∈ R d . Let consider the following quantity for each Borel Λ

L Λ = ∑ ξ ∈X∩Λ f (ξ ,t + ξ , X \ {ξ }).
It is seen that L Λ is covariante:

L τ y Λ (τ y x) = L Λ (x)
and additive: if

Λ 1 ∩ Λ 2 = / 0, L Λ 1 ∩Λ 2 = L Λ 1 + L Λ 2 .
Hence, from the Nguyen and Zessin [START_REF] Nguyen | Ergodic theorems for Spatial Process[END_REF] ergodic theorem, and additionally assume that P is ergodic, C n (t) converges almost surely to

L = 1 |Λ 0 | ∑ ξ ∈X∩Λ 0 f (ξ ,t + ξ , X \ {ξ })
where Λ 0 is a bounded Borel subset of R d . It follows immediately from the GNZ formula (2.1)

L = 1 |Λ 0 | E ∑ ξ ∈X∩Λ 0 f (ξ ,t + ξ , X \ {ξ }) = 1 |Λ 0 | E R d
11 Λ 0 (ξ ) f (ξ ,t + ξ , X)λ (ξ , X)dξ .

Since X is stationary with translation invariant interaction function and using the finite range property (3.5), i.e.

λ (ξ , X) = λ (ξ , X ∩ B(ξ , R)) = β when X ∩ B(ξ , R) = / 0 then, it follows that L = β C(o,t). We conclude that C n (t) defined by ( 4) turns out to be unbiased estimator of β C(o,t) and strongly consistent as n tends infinity.

Using the monotony of functions C and C n , we can approach the functions C and C n by their values in a finite number of points, furthermore by the result (4.15), we have as n → ∞,

  3, . . .} and any Λ ∈ B d 0 . ∑ = denotes summation over summands with index tuples having pairwise distinct components. The origin in R d will be denoted by o.

  ). The sequence of bandwidths (b n ) n≥1 is a decreasing sequence of positive real numbers satisfying: lim n→∞ b d n = 0 and lim n→∞ b d n |Λ n | = ∞.

  sup

  t∈K 0 C n (t)β C(o,t) -→ 0 (4.16) almost surely. By Corollary 2, the result (4.16) and the expression (4.14), we complete the proof.

Acknowledgment

The research was supported by laboratory Jean Kuntzmann, Grenoble Alpes University, France.

Proof of Corollary 2. Using the second-order Georgii-Nguyen-Zessin formula (2.2) with

we have

We remember the second order Papangelou conditional intensity by:

Using the finite range property (3.5) for each function λ (ξ , X) and λ (η, X ∪ {ξ }), we have

We obtain by stationarity of X and from the definition of f giving by (3.1) and C(o,t) is defined through (3.6):

In this way we get: