Mapping the intensity function of a non-stationary point process in unobserved areas

Seismic networks provide data that are used as basis both for public safety decisions and for scientific research. Their configuration affects the data completeness, which in turn, critically affects several seismological scientific targets (e.g., earthquake prediction, seismic hazard...). In this context, a key aspect is how to map earthquakes density in seismogenic areas from censored data or even in areas that are not covered by the network. We propose to predict the spatial distribution of earthquakes from the knowledge of presence locations and geological relationships, taking into account any interaction between records. Namely, in a more general setting, we aim to estimate the intensity function of a point process, conditional to its censored realization, as in geostatistics for continuous processes. We define a predictor as the best linear unbiased combination of the observed point pattern. We show that the weight function associated to the predictor is the solution of a Fredholm equation of second kind. Both the kernel and the source term of the Fredholm equation are related to the first-and second-order characteristics of the point process through the intensity and the pair correlation function. Results are presented and illustrated on simulated non-stationary point processes and real data for mapping Greek Hellenic seismicity in a region with unreliable and incomplete records.


Introduction
Mapping is a key issue in environmental science. A common and first example lies in ecology when mapping species distribution. When the location of individuals is known, we estimate the local density (usually by kernel smoothing), the so-called intensity in point process theory. However, point locations are usually unreachable at the survey scale, so that sampling methods are used; distance sampling or quadrat sampling approaches are possibly the most common ones. When no covariate is available, a mean density estimation is then performed.
However, species distribution characteristics vary spatially as they are governed by environmental data. Several approaches have been developed in that way for species data formed by reported presence locations, also called occurrence-only records (pure records of locations where a species occurred). Generally called Species Distribution Models (SDM), they aim to explain species occurrences from environmental variables. If they are used to gain ecological and evolutionary insights, they are also widely used for model-based interpolation across landscapes or to predict distributions to new geographic regions (Elith and Leathwick, 2009). SDMs are often based on multivariate statistical analysis methods, such as Generalized Linear/Additive Models (GLM/GAM). The most popular models are Maxent (Phillips et al., 2006) and Maxlike (Royle et al., 2012). In these models point locations are aggregated to grid cells and whether one or several individuals are observed in a cell, a one is recorded. Then, their aim is to estimate occurrence probability (species' probability of presence in a grid cell) maps. See e.g. Merow and Silander (2014) for a comparison of the two models and recommendations about their use.
Point process models offer a natural framework for species distribution modelling. Key concerns about SDMs lie in the loss of significant information about the spatial distribution when aggregating point locations, and the dependence of the results on the spatial resolution (Renner and Warton, 2013). Although point process models are connected to Maxent (Aarts et al., 2012;Renner and Warton, 2013), they use a continuous landscape rather than a discretized one, and the number of records is observed and comes from a random process rather than fixed (number of cells/quadrats). Renner et al. (2015) showed that using point process models presents many advantages, including some clarification about the response variable and model assumptions which, in addition, can be checked. Furthermore, because they operate at the individual level, point process models can incorporate interaction between individuals and dependence to environmental covariates.
Another concern about all these approaches is that they are Poisson model-based: their intrinsic definition does not account for relationship between individuals. These interactions nevertheless exist. Competition among individuals often leads to empty areas around each point, so-called exclusion by distance, mimicked by inhibition models. The American Redstarts, for example, compete with conspecifics for habitat in their winter grounds: the re-occupation of the empty areas supports the hypothesis that territoriality in this species acts to exclude conspecifics from certain winter habitat (Marra et al., 1993). On the contrary, individuals can be arranged by groups as with gregarious animals, such groups can also describe some local dispersion of the species around parents (as with plants). This arrangement is achieved in cluster models. The Shorea congestiflora is a dominant species in a 25-ha forest dynamics plot in a rain forest at Sinharaja (Sri Lanka), which apparently shows clustering at several scales (Wiegand et al., 2007). These effects can be mixed, with individuals arranged in groups, but at certain distance of each other inside each group. This can be the case for the spatial distribution of Northern Gannet (Chadoeuf et al., 2011). Similar issues arise in environmental science. The case we consider in this paper concerns seismic mapping in the Aegean Sea region. This is an area with high seismic density due to the presence of numerous volcanoes and plate movements. Earthquakes have a heterogeneous spatial distribution and we might be interested in understanding why certain regions are more favorable than others. Modeling earthquakes will therefore require to take into account geological information and interactions between events, clustering being often linked to aftershocks. Evaluating the seismic hazard requires a reliable monitoring network and sufficient coverage, what is rarely the case, and a lack of recording may not be due to the absence of an earthquake, but to an insufficient or unreliable network. The question then arises of mapping seismic activity in areas where the observation is unreliable or even when the area is not covered by the network, in order to access a relevant map where the seismicity is particularly high.
From a statistical point of view, the method developed in Gabriel et al. (2017) aimed to pre-dict the local variations in the intensity of a spatial point process accounting for the individual relationships modeled by the pair correlation function (which is related to the probability to find a second point of the process at a given distance from a known point of the process). The main interest of this approach is that it estimates, using only first-and second-order characteristics of the point process, the local intensity outside the observation window, hereafter called prediction. The prediction of the local intensity is obtained conditionally to the records in the observation window. However, this method did not allow to consider environmental covariates and thus did not take into account potential spatial variations driven by these covariates at large scale, which may lead to unrealistic predictions in ecology and environment.
To fill this potential unwished situation, we propose to predict the spatial distribution of earthquakes from the knowledge of presence locations and geological relationships, taking into account any interaction between records. Namely, in a more general setting, we aim to estimate the intensity function of a point process, conditional to its censored realization, as in geostatistics for continuous processes. We define a predictor as the best linear unbiased combination of the observed point pattern, where the weight function associated to the predictor is the solution of a Fredholm equation of second kind related to the first-and second-order moments of the point process. We describe our approach in Section 2. We evaluate the goodness of our predictions through a simulation study (Section 3) for several cluster models. In Section 4, our methodology is applied to predict and map Hellenic seismicity in a region with unreliable and incomplete records.

Predicting the intensity conditionally to the observation
We consider in the following a spatial point process Φ in R 2 , i.e. a random pattern of points for which both the number of points and their locations are random. Let us denote Φ(B) the number of points of Φ in any Borelian set B, and Φ B their locations in B. We assume that the point process Φ is simple (i.e. the probability of all points of Φ being distinct is one) and that it has a density probability with respect to the unit rate Poisson process Y .
The intensity function, denoted by λ(x), is defined as the function such that B λ(x) dx = E [Φ(B)], for B any Borelian set; this corresponds to the local probability to observe a point of Φ at a fixed location (if dx is an infinitesimal volume around location x, then P [Φ( dx) = 1] / dx = λ(x) as dx → 0). The intensity function provides a trend in the spatial variation of points density, and we suppose that the intensity is driven by spatial covariates Z, such that λ(x) = h(Z x ), where both h and Z are known. The interaction between points is described through the pair correlation function g(x, y), which gives the extent to which the probability to find a point at a location y changes by the presence of a point of the process at location x (if dy is an infinitesimal volume around y, then P [Φ( dy) = 1 | x ∈ Φ] / dy = g(x, y)λ(y) as dy → 0). We assume that the process is second-order intensity reweighted stationary (Baddeley et al., 2000). This assumption means that its intensity varies in space, but the pair correlation function between two locations depends only on their difference vector.
Here we consider W ⊂ R 2 , a window of interest, and we assume that Φ has only been observed in some observation window W obs ⊂ W . To predict the remaining point process Φ ∩ {W \W obs }, Coeurjolly et al. (2017) express the conditional distribution of Φ∩{W \W obs } given Φ∩W obs = Φ W obs in terms of the conditional density w.r.t Y W \W obs and where f W obs is the marginal density of Φ ∩ W obs w.r.t Y W obs . Unfortunately, the density of Φ restricted to W obs is barely tractable (and hard to handle) except for just some few processes, such as Poisson, Gibbs and determinantal processes; however, this is not the case for Cox and cluster point processes that are often used to model environmental or ecological point patterns.
Our aim is not to predict Φ ∩ {W \W obs } but rather to estimate its local intensity variations at any location x o ∈ W \W obs conditionally to h, Z and Φ W obs . Following Gabriel et al. (2016Gabriel et al. ( , 2017, we refer to it as the spatial local intensity of Φ , that we define by the following limit Hereafter, we denote the conditioning Φ ∩ W obs = Φ W obs , by "|Φ W obs ". We predict the local intensity at points x o ∈ W \W obs using a linear predictor mimicking kriging, it can be written as where δ denotes the Dirac delta function and w is a weight function solution of the Fredholm equation of the second kind (see Apprendix A for details) : with kernel and source term The solution of this Fredholm equation is implicit but it can be numerically approximated using the Galerkin approximation method (Kress, 2013). Briefly, let T h be a given mesh partitioning W obs , e.g obtained by triangulation, and consider {ϕ i } i=1,...,N a basis of the approximation space V h , N = dim V h . We approximate the weight function as a combination of the basis elements : . This approximation plugged into the Fredholm equation leads to a linear problem for all ϕ i : Introducing a matrix formulation of the previous equation reformulates as the so-called Galerkin using some projection on the approximation space K = (K lm ) l,m and K = M KM . We finally get the weights from Equation (3). In all the illustrations that follow we use this approach to approximate the weight functions and make the prediction.

Simulation study
The objective of this section is twofold. First, we want to visualize how the prediction is affected by both the distance from the prediction point x o to the observed window W obs and the point process structure, and we want to measure the variability between our estimator λ(x o |Φ W obs ) and the conditional intensity λ(x o |Φ W obs ). Second, we want to analyze the sensitivity of the predictions when both the intensity and the pair correlation function are unknown, as it is the case in practice.

Goodness of prediction
In order to see how the prediction estimator behaves with respect to the pattern structure we use the inhomogeneous Matérn cluster model which allows the computation of its conditional intensity. This model is obtained as follows. We first define a stationary Matérn cluster process Φ M at , i.e. a process where each point of a Poisson parent process is replaced by a Poissonnian cluster of offspring uniformly distributed in a disc of radius R around the parent point. We then focus on the independent p(x)−thinned process, where p(x) is a deterministic function on R 2 with 0 ≤ p(x) ≤ 1.
Every point x belonging to Φ M at is deleted with probability 1 − p(x), and again its deletion is independent of locations and possible deletions of any other points (Chiu et al., 2013). Let Φ be the process of the thinned offspring and Φ p the parent process (Poisson process with intensity κ). The process Φ is second-order intensity reweighted stationary with intensity λ(x) = κµp(x), with µ the expected number of offspring per parent. The pair correlation function is the one of Φ M at : According to Gabriel and Chadoeuf (2021), the local intensity of the thinned Matérn cluster process Φ is where ∂W stands for the outside border of thickness R of the observation window, b(x o , R) the disc of center x o and radius R, and λ p (y|Φ W obs ) is the conditional intensity of parents in W obs ∪ ∂W given the realization of offspring in W obs . This intensity is approximated by where the first term is the empirical intensity of parents given the observed offspring, and the second term is the intensity of parents with unobserved offspring. See Gabriel and Chadoeuf (2021) for the approximation of the conditional distribution of parent points given the offspring points and its validation.
The inhomogeneous Matérn Cluster process (IMCP) Φ depends on four parameters: the thinning probability p(x), the intensity of parents κ, the mean number of points per parent µ, and the radius of dispersion of the offspring around the parent points R. In our simulation study we fix κ = 50 and µ = 40, and we consider: • the unit square as study region W . The observation window is The size and shape of the prediction windows W pred have been designed to emphasize any effect of the thinning probability on the predictions. We perform n = 1000 simulations of Φ for all combinations of pairs (p(x), R). Each scenario is denoted IMCP(p(x), R). Realizations of IMCP(p 1 (x), R) (resp. IMCP(p 2 (x), R)) are given in the first column of Figure 1 (resp. Figure 2) for the different values of R, illustrating inhomogeneous patterns with increasing range of clustering from top to bottom. For each simulation, we compute the local intensity from both the predictor (1) and the approximation (4) in W pred . We consider the same mesh for the Galerkin approximation method for all configurations, for which W obs is subdivided in 15,194 triangles. The local intensity and the prediction are respectively plotted in the second and third columns of Figures 1 and 2. For a visualization purpose, we plotted the local intensity in W pred (central square/rectangle), as well as a Gaussian kernel smoothing of the intensity in W obs . Zoom of the local intensity and the predictions in W pred are given in Figures 11 and 12 (see Appendix B).
We can see that the method reproduces the structures of the point process. In particular, it reproduces clusters, as soon as there are points close enough to the boundary of W pred , and the prediction tends to the intensity of the process, λ(x) = κµp(x), when it is made at distances larger than R from the boundary of W obs . This is particularly clear for small values of R (first rows in Figures 1 and 2). Results from these realizations show, however, some weak differences between the local intensity and our prediction that we now quantify.
We measure the precision and the variability of our prediction through the relative bias (RB) and the relative root mean square error (RRMSE) as follows Figure 1: First column: inhomogeneous Matérn cluster point patterns, Φ W obs , obtained from p 1 (x). Second column: approximation of the local intensity in W pred (central square) and kernel smoothing of the pattern in W obs . Third column: predicted local intensity in W pred (central square) and kernel smoothing of the pattern in W obs .
where λ i and λ i stand for the prediction and the local intensity of the i-th simulation. We compute RB(x o ) and RRMSE(x o ) at all x o ∈ W pred . This allows us to get an insight of their distribution as illustrated in Figure 3a and 3c for all patterns IMCP(p(x), R). We also look over the values of RB(x o ) and RRMSE(x o ) according to the distance d(x o , W pred ) between the point x o where the prediction is made and the boundary of the prediction window W pred , as plotted in Figure 3b and 3d for all patterns IMCP(p(x), R).
In all cases, the relative bias is concentrated around 0 and is much less than 5% (in absolute value) which reveals an excellent precision of the predictor. The distributions of RRMSE(x o ) show a decreasing variability as the interaction radius R increases. Indeed, as shown in Figures 1 and 2, the inter-point distances of the point patterns are greater for larger values of R and the local intensity is thus more diffuse for large values of R, and more picky for very small values. The weight function is closely related to the pair correlation function. It shows a radial wavy behavior around x o with positive decreasing values at distances less that R, and it has negative values between R and 2R, . . . , and then is null (see Figure 4). As the predictor is the sum of the weights at x ∈ Φ W obs , it reflects this behavior and clusters are thus well identified. The approximation of the local intensity (4) is smoother as it is only expressed in terms of first-order moments. Hence, we get more variability in the border ∂W , and even higher for picky intensities as for R = 0.05. That is also shown in Figures 3b and 3d which illustrate the relative bias and relative RMSE computed at different distances from ∂W . They are both higher at very short distances.

Sensitivity to the estimation of the first-and second-order moments
In the previous section the intensity and the pair correlation functions of the inhomogeneous Matérn process were assumed to be known. We now aim to measure the effect of their estimation on the predictions. We focus on the process IMCP(p 1 (x), 0.09) and we consider four situations for the prediction:   • using the theoretical intensity and pair correlation functions (as defined in the previous section). The resulting predictions are used as the basis for further comparisons.
or the exponential model: g exp (r) = 1 + 1 α 3 exp (−α 4 √ r). This model allows us to check the effect of a misspecification of the pair correlation function.
or non-parametric estimation, g emp , obtained by kernel smoothing, as implemented in the pcfinhom function in spatstat with f λ (x) as intensity function. This situation is often the most natural one and also the one with the less hypotheses. Because empirical estimates do not always satisfy theoretical conditions of a pair correlation function (positive definiteness, consistence), we made an a posteriori selection of the admissible fitted pair correlation functions (see Appendix C).
All results are presented for 250 admissible simulations. We consider a maximum likelihood estimator of the intensity of a Poisson point process to estimate parameters β 1 and β 2 , i.e.
Parameters α i , i = 1, . . . , 4 are estimated by non-linear least squares. Fitted pair correlation functions are plotted in Figure 5. Summaries of fitted parameters are postponed to Appendix C (Table 2). To compare the predictions we compute RB(x o ; λ, λ) and RRMSE(x o ;λ, λ), where λ denotes the prediction using either the theoretical moments or the fitted ones and λ denotes the local intensity (4). As in the previous section, the relative bias and relative RSME are computed at all x o ∈ W pred so that we can plot their distribution, see Figure 6a and 6b. These figures show that the relative bias is about one tenth of the relative RMSE. The relative bias tends to be positive, slightly larger when the predictions are made from (f λ , g exp ). The relative RMSE are in the same order of magnitude, slightly smaller for (f λ , g exp ).
To emphasize the effect of using estimated moments rather than the theoretical ones, we also compute RB(x o ; λ, λ theo ) and RRMSE(x o ; λ, λ theo ), where λ (resp. λ theo ) denotes the prediction using the fitted intensity and pair correlation functions (resp. the theoretical ones). Figure 6c and 6d illustrate the boxplot of their values over all x o ∈ W pred . The ratio between the relative bias and relative RMSe remains the same. The relative bias is close to zero for both the predictions made using (f λ , g emp ) and (f λ , g mat ) and slightly larger and positive for (f λ , g exp ). The relative RMSE is smaller for (f λ , g emp ).
In the former case, the relative bias and relative RMSE are computed from the local intensity which is a slightly smoothed version of the true conditional intensity, whereas in the later case they are computed from the prediction obtained with the theoretical intensity and pair correlation functions leading to less smoothed predictions. As the exponential model also smoothes the predictions, its relative RMSE in the former case is low, but increases in the later case and become similar to the one of from the Matérn model of the pair correlation function. Both predictions from (f λ , g emp ) and (f λ , g mat ) have similar behavior than the one from the theoretical moments. Figure 7 compare the predictions from the different cases on a single simulation and illustrate all these comments. Whilst the results from the non-parametric estimation of the pair correlation function are promising, they do not systematically rely on an admissible pair correlation function (leading to numerical instabilities in the predictions).

Predicting the density of earthquakes
In this section we focus on the Greek-Hellenic area, a region of high seismic hazard due to both tectonic and volcanic seismogenic sources. This seismicity is a result of some motion-induced deformations: northward motion of the African lithosphere, westward motion of the Anatolia plate and collision between African and Eurasian plates (Papazachos and Comninakis, 1971;Le Pichôn and Angelier, 1979;McKenzie, 1972;Anderson and Jackson, 1987). A total of 1173 earthquakes of magnitude greater or equal than 4 occurred in the study area W (black square in Figure 8  Unified Seismological Network (HUSN). Seismic networks provide data that are used as the basis both for public safety decisions and for scientific research. They are essential tools for observing earthquakes and assessing seismic hazards that can be described and characterized to assess their degree of coverage (Siino et al., 2020). Indeed, their configuration affects the data completeness, which in turn, critically affects several seismological scientific targets (e.g., earthquake prediction, seismic hazard, etc. (Vamvakaris et al., 2013)). From different indicators (magnitude of completeness, Radius of Equivalent Sphere that estimates an average error of location), D' Alessando et al. (2011) identified some seismogenic areas that probably are not adequately covered by the HUSN. Based on their results we delineate a zone (in light yellow in Figure 8) with unreliable or missing records. We removed the 48 points located in this yellow zone to predict the local intensity of earthquakes.
Modelling earthquakes needs to account for geological information and interactions between earthquakes, usually in terms of clustering. The most popular models of seismological events are the so-called self-exciting, such as Hawkes and ETAS models (Adamopoulos, 1976;Ogata, 1988). Our approach, however, is model-free, it accounts for environmental heterogeneity and interaction between events. Following Siino et al. (2017), we assume that the intensity of earthquakes loglinearly depends on several spatial geological covariates: D pb , the distance to the plate boundary (orange dashed curves in Figure 8); D f , the distance to the nearest fault (blue lines in Figure 8); and D v , the distance to the nearest main volcano (red triangles in Figure 8). Thus, the intensity takes the form (5) with φ 1 = 6.73, φ 2 = 43.48, φ 3 = 54.783 and φ 4 = 112. The estimated values of model parameters are obtained by the method of maximum likelihood from all points W obs and are given in Table 1, and the estimated intensity is depicted in Figure 9 (left panel, log-scaled). This figure shows that the region to be predicted (bordering in black) is across zones of strong and weak intensities. We Table 1: Estimated parameters of the intensity model λ(x).

Parameter
β 0 β 1 β 2 β 3 β 4 β 5 β 6 Estimate -178.483 0.611 9.673 0.034 -0.061 -0.114 -0.034 Parameter β 7 β 8 β 9 β 10 β 11 β 12 Estimate -0.0268 -0.0175 -0.0179 -0.0214 0.0040 -0.0093 then calculated the empirical pair correlation function using kernel methods and the parametric estimate of the intensity (see Figure 9 (right panel)). We then fitted an exponential model of the form g(r) = exp(−α 2 √ r)/α 1 , withα 1 = 8.9502 andα 2 = 0.0266 (as it is shown in Figure 9 (right panel), for the theoretical fitted model). Model parameters are estimated as in section 3.2. The high values of the pair correlation function at short distances are indicative of a cluster process with a small range, less than 50 km. Then we have g(r) = 1 (horizontal dark grey line), indicating there is no particular spatial structure or interaction. Note that g(0) > 7.5 corresponds to a high probability of observing an earthquake within a small disc centered at an observed earthquake.
The prediction of the local intensity is obtained by using a mesh of 22,686 triangles for the Galerkin approximation method. The left panel of Figure 10 shows the predicted local intensity λ(x o |Φ W obs ) in W pred = W \W obs and a Gaussian-kernel smoothing of observed earthquake locations in W obs , with bandwidth 22 km. Both intensities are at a log-scale. The blue area in the center of W obs indicates a kernel smoothing value near zero because of the absence of points in this region (see Figure 8). Compared with Figure 9, this plot also emphasizes the differences between the intensity fitted from covariates in Equation (5) and the empirical intensity obtained by kernel smoothing. This leads to conclude that we can expect some influence from point locations on the local intensity. This effect is observed in the middle panel of Figure 10. It focuses on the prediction window in which both the predicted local intensity and the fitted intensity of the point pattern λ(x o ) are plotted at a log-scale. The right panel plots their ratio to highlight their differences. The different range values between λ(x o |Φ W obs ) and λ(x o ) indicate that if we did not take into account the effect of earthquakes locations in W obs , we would get regular variations of the local intensity in W pred . The knowledge of these points add local hot spots close to the border between W obs and W pred . Hence, the prediction window shows a high and spatially varying local intensity, whose variations reflect the scales of the structures of the underlying process, through both the intensity and the pair correlation function. Note that from a practical point of view, the approximation of the weight functions has been implemented in FreeFem++ (Hecht, 2012) a high level integrated development environment for numerically solving partial differential equations. All other computations are implemented in R (R Core Team, 2021). We used the R package spatstat (Baddeley and Turner, 2005) for fitting the intensity λ(x) and for computing the empirical pair correlation function. Because interfacing R and FreeFem++ is not straightforward and not feasible on any operating system, codes are only available upon request to the first author.

Discussion
We have proposed a predictor of the local intensity of point processes conditionally to some censored observations, and accounting for the individual relationships and considering environmental covariates. In this general context we can distinguish two situations: in a first case, censoring is a continuous function defined over the study region W , so that observed and unobserved points share the same space; in a second case, observations are taken in some windows and we want to predict outside the observation window. In this paper we presented the latter case and considered that the spatial variation driving the point density is related to an intensity governed by environmental covariates. When observed and unobserved points share the same space, the study region is the observation window and the observation is not exhaustive. In the example of earthquakes, we know the network reliability. In some sense it can be viewed as a probability of detection, say π(x), which is independent of the underlying process of earthquakes. This process is the union of the observed process, Φ π , and the unobserved process, Φ 1−π . Then the local intensity is the sum of the intensity of the observed process and the intensity of the unobserved process given the observation, what can be predicted as previously. But now, weights also depend on the observation probability.
The interaction with other processes provides another form of spatial variation. We can easily imagine that the presence or the absence of a species depend on the presence/absence of other species. The extension of our approach to multi-type processes is an on-going work. In that case, we consider several processes, say Φ (k) , k = 1, . . . , K, observed respectively in some windows W 1 , ..., W K . To predict the intensity of the first process given the realization of the others, we can define the predictor as follows The unbiasedness constraint on the predictor modifies the constraint on the weight function, and minimizing the error prediction variance under this constraint leads to a new Fredholm equation which now depends on the pair correlation function of each process and on the cross-pair correlation functions.
A question remains: how to control the approximation accuracy when solving the Fredholm equation? There are two ways to increase the accuracy of the computed solution of the Fredholm equation: either refine the mesh, i.e. lower the size of the triangles, or use another approximation basis than the usual one (P1={piece-wise Linear functions}). Both options significantly increase the dimension of the resulting dense linear problem and from the point of view of computing efficiency we can wonder which is the most suitable approximation. Note that for a given mesh, we get a more precise representation of the regularity of the solution if we use a finite element basis involving more regular functions e.g. P2={piece-wise Quadratic functions} or P3={piece-wise Cubic functions}. . . However, changing the basis of approximation won't reduce the approximation errors made with respect to the true solution of the Fredholm equation, errors that are intrinsically linked to the typical size of a triangle of the mesh. Up to our knowledge, no rule has been published that relates the size of the mesh to the characteristics of the equation. A too crude size can even lead to negative estimated local intensities. We propose to verify empirically the approximation quality by simulating point realizations using an easy to manipulate point process with the same first-and second-order characteristics, as the log-Gaussian Cox point process, and check discrepancy between reestimated local intensities and the estimated one.
The constraint over the spatial weights in order to obtain an unbiased prediction can be expressed as then under non-stationary assumption it becomes On the other hand, the predictor must minimize the error prediction variance, This is equivalent to minimise the following expression We consider the Lagrange multipliers under the constraint for the spatial weight function and we define Finally, we can rewrite the previous expression as and therefore Using variational calculation and the Riesz representation theorem, it follows that and therefore Considering the integral of previous equation over the spatial window W obs and respect to x, it follows that and then we obtain from which we get Finally, plugging µ/2 into Equation (6) λ(x o |Φ w obs ) λ(x o |Φ w obs ) Figure 11: Conditional intensity (left panels) and prediction of the local intensity (right panels) obtained from IMCP(p 1 (x), R) with R = 0.05 (first row), R = 0.09 (second row) and R = 0.013 (third row).

C First-and second-order moments estimation
In section 3 we analyze the effect of plugging estimates of the intensity and the pair correlation functions in the predictor. We considered a parametric model for the intensity function, f λ (x), two parametric (Matérn and Exponential) and one non parametric (referred to as empirical and denoted g emp (r)) models for the pair correlation function. Table 2 summarizes the fitted parameters of the different models. The non parametric estimate of the pair correlation function is obtained by kernel smoothing, as implemented in the pcfinhom function in spatstat with f λ (x) as intensity function. However this function may lead to odd estimates that often do not satisfy theoretical conditions of a pair correlation function (positive definiteness, consistence). Odd estimates are plotted in Figure 13 and show that they do not tend to 1 as r increases. Regarding the expression of the weights (Equation 2), this highly impact the predictions: we get very hight values for the related simulations. For instance from the estimates plotted in Figure 13 we get predictions between −10130 and 13295. Hence for the simulation study, we made an a posteriori selection of the admissible fitted pair correlation functions as follows : • median {|g emp (r) − 1|;r > 2R} < 0.05, i.e. the median distance between g emp (r) and 1 for all r > 2R must be less than 0.05, • the predictions must be positive and not exceed 6000.