Jean-François Coeurjolly
email: coeurjolly.jean-francois@uqam.ca

Frédéric Lavancier
email: frederic.lavancier@univ-nantes.fr

Approximation intensity for pairwise interaction Gibbs point processes using determinantal point processes

Keywords: MSC 2010 subject classifications: Primary: 60G55, secondary: 82B21 Determinantal point process, Georgii-Nguyen-Zessin formula, Gibbs point process, Laplace transform

The intensity of a Gibbs point process is usually an intractable function of the model parameters. For repulsive pairwise interaction point processes, this intensity can be expressed as the Laplace transform of some particular function. Baddeley and Nair (2012) developped the Poissonsaddlepoint approximation which consists, for basic models, in calculating this Laplace transform with respect to a homogeneous Poisson point process. In this paper, we develop an approximation which consists in calculating the same Laplace transform with respect to a specific determinantal point process. This new approximation is efficiently implemented and turns out to be more accurate than the Poisson-saddlepoint approximation, as demonstrated by some numerical examples.

Introduction

Due to their simple interpretation, Gibbs point processes and in particular pairwise interaction point processes play a central role in the analysis of spatial point patterns (see [START_REF] Van Lieshout | Markov Point Processes and Their Applications[END_REF]; [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]; [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]). In a nutshell, such models (in the homogeneous case) are defined in a bounded domain by a density with respect to the unit rate Poisson point process which takes the form

f (x) ∝ β |x| u∈x g(v -u),
where x is a finite configuration of points, where β > 0 represents the activity parameter, |x| is the number of elements of x and where g : R d → R + is the pairwise interaction function.

However, many important theoretical properties of these models are in general intractable, like for instance the simplest one, the intensity λ ∈ R + , representing the mean number of points per unit volume. It is known (see e.g. Section 2.2) that

λ = β E u∈x g(u) .
Such an expectation is in general intractable. As clearly outlined by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF], this intractability constitutes a severe drawback. For example, simulating a Gibbs point process with a prescribed value of λ cannot be done beforehand even for simple models such as Strauss models. [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF] suggest to evaluate the expectation with respect to a homogeneous Poisson point process with intensity λ. This results in the Poisson-saddlepoint approximation, denoted by λ ps , obtained as the solution of log λ ps = log β -λ ps G where G = R d (1 -g(u)) du (provided this integral is finite). The general idea of the present paper is to evaluate the same expectaction with respect to a determinantal point process (with intensity λ). Determinantal point processes (DPP), see e.g. [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], are a class of repulsive models which is more tractable than Gibbs models. For example all moments are explicit. If g ≤ 1 and has a finite range R > 0, our approximation denoted by λ dpp is the solution of

log λ dpp = log β + (1 + λ dpp G/κ) log 1 - λ dpp G 1 + λ dpp G/κ .
where

κ = max |B(0, δ)| (1 -g) 2 ,
(1 -g) 2 |B(0, R)| , imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017

Coeurjolly and Lavancier/Intensity approximation using DPP

3

|A| denotes the volume of some bounded domain A ⊂ R d , B(0, ρ) is the Euclidean ball centered at 0 with radius ρ and δ ≥ 0 is some possible hard-core distance.

Both approximations λ dpp and λ ps can be obtained very quickly with a unitroot search algorithm. Figure 1 reports λ dpp and λ ps as well as the true intensity λ (obtained by Monte-Carlo simulations) for Strauss models in terms of the interaction parameter γ 1 ∈ [0, 1]. This setting is considered by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF]. The DPP approximation outperforms the Poisson-saddlepoint approximation especially when γ 1 is close to zero, i.e. for very repulsive point processes. More numerical illustrations are displayed in Section 4. The rest of the paper is organized as follows. Section 2 provides necessary notation and background material on point processes, Gibbs point processes and determinantal point processes. Intensity approximations are discussed in detail in Section 3. Finally, Section 4 presents numerical experiments for several classes of pairwise interaction point processes. (2003) and [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF].

Gibbs point processes and determinantal point processes

Background and Poisson point processes

A spatial point process is said to have an nth order intensity function ρ (n) if for any nonnegative measurable function h : (R d ) n → R + , the following formula referred to as Campbell-Mecke formula holds

E = u1,...,un∈X h(u 1 , . . . , u n ) = R d . . . R d h(u 1 , . . . , u n )ρ (n) (u 1 , . . . , u n ) du 1 . . . du n , (2.1)
where the sign = over the sum means that u 1 , . . . , u n are pairwise distinct. Then, ρ (n) (u 1 , . . . , u n ) du 1 • • • du n can be interpreted as the approximate probability for X having a point in each of infinitesimally small regions around u 1 , . . . , u n of volumes du 1 , . . . du n , respectively. We also write ρ(u) for the intensity function ρ (1) (u). A spatial point process X in R d is said to be stationary (respectively isotropic) if its distribution is invariant under translations (respectively under rotations). When X is stationary, the intensity function reduces to a constant denoted by λ in the rest of this paper. As a matter of fact, λ measures the mean number of points per unit volume.

The Poisson point process, often defined as follows (see e.g. [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]), serves as the reference model. Definition 2.1. Let ρ be a locally integrable function on S, for S ⊆ R d . A point process X satisfying the following statements is called the Poisson point process on S with intensity function ρ:

• for any m ≥ 1, and for any disjoint and bounded B 1 , . . . , B m ⊂ S, the random variables X B1 , . . . , X Bm are independent; • N (B) follows a Poisson distribution with parameter B ρ(u) du for any bounded B ⊂ S.

Among the many properties of Poisson point processes, it is to be noticed that the nth order intensity function writes ρ

(n) (u 1 , . . . , u n ) = n i=1 ρ(u i ), for any pairwise distinct u 1 , . . . , u n ∈ S.
Let Z be a unit rate Poisson point process on S, which means that its intensity is constant and equal to one. Assume, first, that S is bounded (|S| < ∞). We say that a spatial point process X has a density f if the distribution of X is absolutely continuous with respect to the one of Z and with density f . Thus, for any nonnegative measurable function h defined on N , E h(X) = E(f (Z)h(Z)). Now, suppose that f is hereditary, i.e., for any pairwise distinct u 0 , u 1 , . . . , u n ∈ S, f ({u 1 , . . . , u n }) > 0 whenever f ({u 0 , u 1 , . . . , u n }) > 0. We can then define imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017 the so-called Papangelou conditional intensity by

λ(u, x) = f (x ∪ u)/f (x) (2.2)
for any u ∈ S and x ∈ N , setting 0/0 = 0. By the interpretation of f , λ(u, x) du can be considered as the conditional probability of observing one event in a small ball, say B, centered at u with volume du, given that X outside B agrees with x. When f is hereditary, there is a one-to-one correspondence between f and λ.

Because the notion of density for Z when S = R d makes no sense, the Papangelou conditional intensity cannot be defined through a ratio of densities in R d . But it still makes sense as the Papangelou conditional intensity can actually be defined at the Radon-Nykodym derivative of P ! u the reduced Palm distribution of X with respect to P, the distribution of X (see [START_REF] Daley | An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods[END_REF]). We do not want to enter in too much detail here and prefer to refer the interested reader to [START_REF] Coeurjolly | A tutorial on Palm distribution for spatial point processes[END_REF].

Finally, we mention the celebrated Georgii-Nguyen-Zessin formula (see [START_REF] Georgii | Canonical and grand canonical Gibbs states for continuum systems[END_REF][START_REF] Nguyen | Ergodic theorems for spatial processes[END_REF], which states that for any h : R d ×N → R (such that the following expectations are finite)

E u∈X h(u, X \ u) = R d E (h(u, X)λ(u, X)) du.
(2.3) By identification of (2.1) and ( 2.3), we see a link between the intensity function of a point process and the Papangelou conditional intensity: for any

u ∈ R d ρ(u) = E (λ(u, X)) ,
which in the stationary case reduces to λ = E (λ(0, X)) .

(2.4)

Gibbs point processes

For a recent and detailed presentation, we refer to [START_REF] Dereudre | Introduction to the theory of Gibbs point processes[END_REF]. Gibbs processes are characterized by an energy function H (or Hamiltonian) that maps any finite point configuration to R ∪ {∞}. Specifically, if |S| < ∞, a Gibbs point process on S associated to H and with activity β > 0 admits the following density with respect to the unit rate Poisson process:

f (x) ∝ β |x| e -H(x) , (2.5)
where ∝ means "proportional to". This definition makes sense under some regularity conditions on H, typically non degeneracy (H(∅) < ∞) and stability (there exists A ∈ R such that H(x) ≥ A|x| for any x ∈ N ). Consequently, configurations x having a small energy H(x) are more likely to be generated by a Gibbs point process than by a Poisson point process, and conversely for configurations having a high energy. In the extreme case where H(x) = ∞, then x cannot, almost surely, be the realization of a Gibbs point process associated to H.

In this paper, we focus on pairwise interaction point processes. To be close to the original paper by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF] the present contribution is based on, we use their notation: a Gibbs point process in S is said to be a pairwise interaction point process with pairwise interaction function g :

R d → R + , if its density writes f (x) ∝ β |x| u,v∈x g(u -v).
If |S| = ∞, this definition and more generally Definition (2.5) do not make sense since H(x) can be infinite or even undefined if |x| = ∞. In this case, Gibbs point processes have to be defined via their conditional specifications and for pairwise interactions Gibbs point processes, restrictions on g have to be imposed for existence (see again [START_REF] Dereudre | Introduction to the theory of Gibbs point processes[END_REF] and the references therein for details). Nonetheless, as mentioned in the previous section, the concept of Papangelou conditional intensity applies whenever |S| < ∞ or |S| = ∞, and in either case it has the explicit form

λ(u, x) = β v∈x g(u -v), (2.6)
for any u ∈ S. Note that when S = R d , a pairwise interaction Gibbs point process is stationary if g is symmetric and it is further isotropic if g(u -v) depends simply on u -v . From (2.4), we deduce that the intensity parameter of a stationary pairwise interaction process writes

λ = E (λ(0, X)) = β E v∈x g(v) .
(2.7)

Let us give a few examples (which are in particular well-defined in R d ). Many other examples can be found e.g. in the recent monograph by [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF].

• Strauss model: let γ ∈ [0, 1] and 0 < R < ∞ g(u) = γ if u ≤ R 1 otherwise. (2.8) • Strauss Hard-core model: let γ ∈ R + and 0 < δ < R < ∞ g(u) =    0 if u < δ γ if δ ≤ u ≤ R 1 otherwise.
• Piecewise Strauss Hard-core model:

g(u) =    0 if u < δ γ i if R i ≤ u ≤ R i+1 , i = 1, . . . , I 1 otherwise,
where

I ≥ 1, 0 ≤ R 1 = δ < R 2 < • • • < R I+1 = R < ∞ and γ 1 , . . . , γ I ∈ R + if δ > 0, otherwise γ 1 , . . . , γ I ∈ [0, 1]. • Diggle-Graton model: let γ ∈ [0, 1] g(u) = u R 1/γ if u ≤ R 1 otherwise,
where for t ∈ (0, 1), t ∞ = 0 and 1 ∞ = 1 by convention.

Let us note that a Strauss model with γ = 0 and radius R is actually a hard-core model with radius R. The Diggle-Graton potential can be found in [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF] in a slightly different parameterization. The one chosen here makes comparisons with the Strauss model easier. For instance, when γ = 0 the model reduces to a Strauss model with γ = 0 and radius R. When γ = 1, the function g grows linearly from 0 to 1. Figure 2 depicts the form of some of the pairwise interaction functions presented above. A Gibbs point process has a finite range R if for any u ∈ R d and x ∈ N , λ(u, x) = λ(u, x ∩ B(u, R)). For pairwise interaction point processes, this property translates to g(u) = 1 for any u ∈ R d such that u > R. All previous models have a finite range R < ∞. An example of infinite range pairwise interaction point process which will not be considered in this paper is the Lennard-Jones model (see e.g. [START_REF] Ruelle | Statistical Mechanics: Rigorous Results[END_REF]; [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]).

Determinantal point processes

Determinantal point processes (DPPs) are models for inhibitive point patterns. We refer to [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF] for their main statistical properties. They are defined through a kernel function K which is a function from S × S to C. A point process is a DPP on S with kernel K, denoted by DPP(K), if for any n, its nth order intensity function takes the form

ρ (n) (u 1 , . . . , u n ) = det[K](u 1 , . . . , u n ),
(2.9)

for every (u 1 , . . . , u n ) ∈ S n , where [K](u 1 , . . . , u n ) denotes the matrix with entries K(u i , u j ), 1 ≤ i, j ≤ n. In particular, the intensity function of DPP(K) is K(u, u).

Conditions on the kernel K are required to ensure the existence of DPP(K). For our purpose, we will only consider DPPs on a compact set. So let us assume that S is compact and suppose that K is a continuous real-valued covariance function on S × S. In this setting, by the Mercer theorem (see [START_REF] Riesz | Functional Analysis[END_REF]), K admits the spectral expansion

K(u, v) = ∞ i=1 λ i φ i (u)φ i (v), ∀u, v ∈ S, (2.10)
where {φ i } i≥1 is an orthonormal basis of L 2 (S) and where λ i , i ≥ 1, are referred to as the eigenvalues of K. Under the above assumptions, DPP(K) exists if and only if λ i ≤ 1 for all i. Due to their tractability, DPPs have many interesting properties. Many of them have been obtained by [START_REF] Shirai | Random point fields associated with certain Fredholm determinants. I. Fermion, Poisson and boson point processes[END_REF], from which we derive the following key-equation used by our intensity approximation.

Proposition 2.2. Let X be a DPP on a compact set S with kernel K. Assume that K is a continuous real-valued covariance function on S × S whose all eigenvalues are not greater than 1. For any function g :

S → [0, 1] E v∈X g(v) = i≥1 (1 -λi ) (2.11)
where λi , for i ≥ 1, are the eigenvalues of the kernel K : S × S → R given by

K(u, v) = 1 -g(u)K(u, v) 1 -g(v).
Proof. Note that

E v∈X g(v) = L X (-log g)
where L X denotes the Laplace transform of X. From Theorem 1.2 in [START_REF] Shirai | Random point fields associated with certain Fredholm determinants. I. Fermion, Poisson and boson point processes[END_REF], for any nonnegative measurable function f on S L X (f ) = Det(I -K)

imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017

Coeurjolly and Lavancier/Intensity approximation using DPP

where Det denotes the Fredholm determinant of an operator and K is the integral operator associated to the kernel

K(u, v) = 1 -exp(-f (u))K(u, v) 1 -exp(-f (v)).
On the other hand, see for instance (2.10) in [START_REF] Shirai | Random point fields associated with certain Fredholm determinants. I. Fermion, Poisson and boson point processes[END_REF],

Det(I -K) = exp - ∞ n=1 1 n Tr ( Kn ) ,
where Tr denotes the trace operator. The result follows from the fact that for any n ≥ 1

Tr ( Kn ) = S n K(u 1 , u 2 ) • • • K(u n , u 1 )du 1 • • • du n = i≥1 λn i .

Intensity approximation

Poisson-saddlepoint approximation

We remind that the intensity parameter of a Gibbs point process, and in particular a pairwise interaction point process satisfies (2.7). The expectaction in (2.7) is to be regarded with respect to P the distribution of the Gibbs point process X. [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF] suggest to replace P by a simpler distribution, say Q, for which the right-hand-side of (2.7) becomes tractable. The Poissonsaddlepoint approximation consists in choosing Π(λ), the Poisson distribution with parameter λ, as distribution Q. As a result, the Poisson-saddlepoint approximation consists in resolving the equation

λ = β E Π(λ) v∈Y g(v) = β E Π(λ) exp v∈Y log g(v) , (3.1)
with the convention that log 0 = -∞ and where, to avoid any ambiguity, we denote by Y a Poisson point process with intensity λ defined on R d and stress also this by indexing the E with the distribution Π(λ). It turns out that if g(u) ∈ [0, 1] for any u ∈ R d , the right-hand side of (3.1) is the Laplace transform of some Poisson functional and equals β exp(-λG) where G = R d (1-g(u)) du, see e.g. Møller and Waagepetersen (2004, Proposition 3.3). As noticed in [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF], this formula extends to more general functions g, provided G > -∞. Hence, the Poisson-saddlepoint approximation, denoted by λ PS in this paper, is defined as the solution of

λ ps = β exp(-λ ps G) ⇐⇒ λ ps = W (βG) G (3.2)
imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017

Coeurjolly and Lavancier/Intensity approximation using DPP

where W is the inverse function of x → x exp(x).

For stationary pairwise Gibbs models with finite range R, and such that λ(u, x) ≤ β (or equivalently such that g ≤ 1), then 0 ≤ G ≤ |B(0, R)|. In this case, [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF] prove, among other properties, that λ ps exists uniquely and is an increasing function of β. From a numerical point of view, λ PS can be very efficiently and quickly estimated using root-finding algorithms.

DPP approximation

Following the same idea as the Poisson-saddlepoint approximation, for a repulsive stationary pairwise interaction point process with pairwise interaction function g ≤ 1 having a finite range R, we suggest to substitute the measure P involved in the expectation (2.7) by the measure Q corresponding to a DPP defined on B(0, R) with some kernel K (to be chosen) and intensity λ, i.e. K(u, u) = λ. Similarly to the previous section, by letting DPP(K; λ) denote the distribution of such a DPP and Y ∼ DPP(K; λ), the DPP approximation of the intensity λ is the solution of

λ = β E DPP(K;λ) v∈Y g(v) .
(3.3) From Proposition 2.2 and in particular from (2.11), this yields the estimating equation

log λ = log β + i≥1 log(1 -λi ),
where the eigenvalues λi of K are related to λ by the relation

K(u, v) = 1 -g(u)K(u, v) 1 -g(v) with K(u, u) = λ.
To complete this approximation, the eigenvalues λi need to be specified.

In the following we choose the eigenvalues λi to be zero except a finite number N of them that are all equal. Given that

i≥1 λi = R d K(u, u)du = R d (1 -g(u))K(u, u)du = λG,
this means that for some N ≥ λG, λi = λG N , for i = 1, . . . , N (3.4) and λi = 0 for i ≥ N + 1. With this choice, the integer N remains the single parameter to choose in our approximation. Note that N ≥ λG is a necessary condition to ensure λi ≤ 1 and so the existence of a DPP with kernel K, but it is in general not sufficient to ensure the existence of the relation between K and K where K defines a DPP. This will be clearly illustrated below when g is the Strauss interaction function. For the choice (3.4), the DPP approximation of the intensity, denoted by λ dpp , becomes the solution of

log λ dpp = log β + N log 1 - λ dpp G N ⇐⇒ λ dpp = β 1 - λ dpp G N N . (3.5)
To motivate (3.4) and how we should set N , assume for a moment that g is the interaction function of a Strauss model with range R and interaction parameter γ ∈ [0, 1], see (2.8). In this case K(u, v) = (1 -γ)K(u, v) for any u, v ∈ B(0, R) and the eigenvalues λ i of K satisfy λi = (1-γ)λ i . In the approximation (3.3), we start by choosing a kernel K with a finite number of non-vanishing eigenvalues λ i that are all equal. In view of λ i = K(u, u)du = λb, where b denotes the volume of B(0, R), this leads to λ i = λb/N for i = 1, . . . , N and N ≥ λb. Note that the latter inequality is necessary to ensure the existence of DPP(K). Going back to λi , this means that (3.4) follows with the necessary and sufficient condition N ≥ λb = λG/(1 -γ) which is greater than λG.

In order to set N precisely for the Strauss model, remember that a homogeneous DPP is more repulsive when its eigenvalues are close to 1, see [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF]; [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF], and at the opposite a DPP is close to a Poisson point process when its eigenvalues are all close to 0. This suggests that in order to make the approximation (3.3) efficient, we should choose λ i close to 1 when the Gibbs process we want to approximate is very repulsive, that is when γ is close to 0. Moreover the eigenvalues should decrease to 0 when γ increases to 1. If λ i = λb/N , this is equivalent to choosing N an integer that increases from λb to infinity when γ increases from 0 to 1. A natural option is thus to choose N as the smallest integer larger than λb/(1 -γ). Our final choice for the Strauss model is therefore N = λb/(1 -γ) , where . denotes the ceiling function, which we may write, for later purposes, N = λG/(1 -γ) 2 .

However, with the latter choice, the function in the right-hand side of equation (3.5) is not continuous in λ, which may lead to none or several solutions to this equation. As a last step in our approximation, we therefore consider the upper convex envelope of this function, ensuring a unique solution to (3.5). This finally leads for the Strauss interaction process to the approximation λ dpp defined as the solution of

log λ dpp = log β + (1 + λ dpp G/(1 -γ) 2 ) log 1 - λ dpp G 1 + λ dpp G/(1 -γ) 2 .
Let us now discuss the case of a general pairwise interaction function g. In this setting, it is in general not possible to relate the eigenvalues λ i of K with the eigenvalues λi of K. Motivated by the Strauss case, we choose λi as in (3.4) where N = λG/κ and κ ∈ [0, 1] is a parameter that takes into account the repulsiveness encoded in g. In general κ must be close to 0 when g is close to 1 (the Poisson case), and close to 1 when g is close to a pure hard-core interaction. We decide to quantify the repulsiveness of the model by b -1 (1 -g) 2 , in agreement with our choice for the Strauss model for which κ = (1 -γ) 2 . Note that for a pairwise interaction g with range R and involving a possible hard-core distance δ, we have |B(0, δ)| ≤ (1 -g) 2 ≤ |B(0, R)| where the left or right equality occurs for a pure hard-core interaction (if δ > 0), a situation where κ must be 1. This leads us to the choice

κ = max |B(0, δ)| (1 -g) 2 , (1 -g) 2 |B(0, R)| . (3.6)
Plugging N = λG/κ into (3.5) and considering the upper convex envelope to ensure the existence of a unique solution, we finally end up with our general DPP approximation being the solution of

log λ dpp = log β + (1 + λ dpp G/κ) log 1 - λ dpp G 1 + λG/κ (3.7) ⇐⇒ λ dpp = W κ (βG/κ) G/κ , (3.8)
where κ is given by (3.6) and where W κ is the inverse function of

x → x 1 -κx 1+x -1-x
. In view of (3.7)-(3.8) and similarly to λ ps , the approximation λ dpp can be very efficiently implemented using root-finding algorithms. We further have the following properties.

Theorem 3.1. Consider a stationary pairwise interaction process in R d with Papangelou conditional intensity given by (2.6) which is purely inhibitory, i.e. g(u) ≤ 1 for all u ∈ R d and with finite range R. Then, λ dpp exists uniquely, is an increasing function of β and is such that λ dpp ≤ λ ps .

Proof. Let f ps and f dpp denote the two real-valued functions given by

f ps (λ) = β exp(-λG) and f dpp (λ) = β 1 - λG 1 + λG/κ 1+λG/κ .
The approximations λ ps and λ dpp are defined by the fixed point equations λ ps = f ps (λ ps ) and λ dpp = f dpp (λ dpp ). Since for any x ∈ [0, 1), log(1 -x) ≤ -x, we have for any λ 0 ≤ f dpp (λ) ≤ β exp(-λG) = f ps (λ).

(3.9)

In particular 0 ≤ lim λ→∞ f dpp (λ) ≤ lim λ→∞ f ps (λ) = 0. In addition, f dpp (0) = β and it can be verified that f dpp is a decreasing function. Hence the solution to (3.7) exists uniquely. The function W κ can also be shown to be increasing on R + for any κ ∈ [0, 1], so we deduce from (3.8) that λ dpp is an increasing function of β. Finally, (3.9) shows that λ dpp ≤ λ ps .

imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017

Numerical study

In this section, we focus on the planar case to investigate the performances of the DPP approximation and compare it with the initial one proposed by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF]. All computations were performed in the R language (R development core team, 2011). The Poisson-saddlepoint approximation as well as the DPP approximation are implemented using root-finding algorithms and in particular we use the R function uniroot for this task.

We have considered 14 different numerical experiments involving Strauss models (S), Strauss hard-core models (SHC), Diggle-Graton models (DG), piecewise Strauss models (PS) and piecewise Strauss hard-core (PSHC) models. The pairwise interaction functions of these models are detailed in Section 2.2. To sum up here are the parameters, that include a continuously varying parameter γ 1 ∈ [0, 1]: For each numerical experiment, we therefore obtain curves of intensity approximation in terms of γ 1 . For γ 1 varying from 0 to 1 by step of 0.05 (the value 0 is exluded for DG models to save time), the true intensity λ is estimated by Monte-Carlo methods. For each set of parameters m realizations of the model are generated on the square [-2R, 1 + 2R] 2 and then clipped to the unit square. That strategy is detailed and justified by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF]. Specifically, the number of points in each realization is averaged to obtain the estimated intensity and its standard error. The simulation results for the Strauss models with β = 50 or 100 were obtained by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF], where m = 10000 realizations were generated and the exact simulation algorithm was used, implemented in the R function rStrauss of the spatstat package (see [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]). For the Strauss models with β = 200, SHC models, PS and PSHC models, we generate m = 1000 replications and use the rmh function in the spatstat package which implements a Metropolis-Hastings algorithm. Even if we use 10 6 iterations of the algorithm, the results may be slightly biased. For the DG models, the R package spatstat provides an exact simulation algorithm (function rDiggleGraton) and for such models we generate 10000 replications when β = 200 and R = 0.025, 0.05 and when β = 50 and R = 0.15. We used 1000 replications when β = 200 and R = 0.075 to save time.

•
All results can be found in Figures 3, 4 and5. Plots provide the same information: we depict intensity approximation λ based on different methods in terms of γ 1 . The dashed curve represents the Poisson-saddlepoint approximation proposed by [START_REF] Baddeley | Fast approximation of the intensity of Gibbs point processes[END_REF] and detailed in Section 3.1. The solid curve is the DPP approximation we propose in this paper and is given by (3.8).

Let us first comment Figure 3 dealing with Strauss models. As expected the Poisson-saddlepoint approximation is not efficient when γ 1 is small, i.e. for very repulsive models. This is very significant in particular for the Strauss hard-core model, see Figure 3 (f). The DPP aproximation we propose is more likely able to capture the repulsiveness of the Strauss models. Figure 4 also clearly shows that our approximation is particulalry efficient and outperforms unamibigously the Poisson-saddlepoint approximation. Note that replications for the Diggle-Graton models are generated using an exact algorithm; so the numerical results seem to be exact, except the slight bias induced by clipping the pattern from [-2R, 1 + 2R] 2 to the unit square.

We finally comment Figure 5. When γ 2 = 0.5, i.e. Figures 5 (a)-(b), the results are very satisfactory. Our approximation is able to approximate λ very efficiently for any value of γ 1 . For Figures 5 (c)-(d), γ 2 = 0 which means that points within a distance comprised between 0.05 and 0.1 are forbidden. Such a parameterization tends to create repulsive clusters. When γ 1 = 1 and δ = 0, such a piecewise Strauss model was called annulus model by [START_REF] Stucki | Bounds for the probability generating functional of a Gibbs point process[END_REF]. This model demonstrates the limitations of our approximation even if when γ 1 is close to zero which means that the model is close to a hard-core process with radius 0.1 our approximation remains satisfactory.

Fig 1 .

 1 Fig 1. Comparison of the exact intensity (small boxplots), the Poisson-saddlepoint approximation (dashed line) and the DPP approximation (solid line) for homogeneous Strauss models with activity parameter β and range of interaction R. Curves and boxplots are reported in terms of the interaction parameter γ 1 ∈ [0, 1].

Fig 2 .

 2 Fig 2. Examples of pairwise interaction functions for the Strauss model (γ = 0.5, R = 0.15), the piecewise Strauss hard-core model (δ = R 1 = 0.05, R 2 = 0.1, R 3 = R = 0.15, γ 1 = 0.8, γ 2 = 0.2), and the Diggle-Graton model (γ = 0.05, 0.3 and 1, R = 0.15).

  Strauss (S): β = 100 with R = 0.05 or 0.1; β = 50 with R = 0.1 or 0.15; β = 200 with R = 0.05. For all these models γ = γ 1 . • Strauss hard-core (SHC): β = 200, δ = 0.025, R = 0.05. For this model γ = γ 1 . • Diggle-Graton (DG): β = 200, R = 0.025, 0.05 or 0.075 and β = 50 and R = 0.15. For all these models γ = γ 1 . • Piecewise Strauss and Strauss hard-core (PS and PSHC): β = 200, δ = 0 or 0.025, γ = (γ 1 , γ 2 ) with γ 2 = 0 or 0.5. The vector of breaks is R = (0.05, 0.1).

Fig 3 .Fig 4 .Fig 5 .

 345 Fig 3. Comparison of the exact intensity (small boxplots obtained by Monte-Carlo method), the Poisson-saddlepoint approximation (dashed line) and the DPP approximation (solid line) for homogeneous Strauss and Strauss hard-core models with activity parameter β, range of interaction R and eventually hard-core distance δ. Curves and boxplots are reported in terms of the interaction parameter γ 1 ∈ [0, 1].

  For d ≥ 1, let X be a spatial point process defined on R d , which we see as a random locally finite subset of R d . Local finiteness of X means that X B = X∩B is finite almost surely (a.s.), that is the number of points N (B) of X B is finite a.s., whenever B ⊂ R d is bounded. We let N stand for the state space consisting of the locally finite subsets (or point configurations) of R d . Let B(R d ) denote the class of bounded Borel sets in R d . For any B ∈ B(R d ), we denote by |B| its Lebesgue measure. A realization of X B is of the form x = {x 1 , . . . , x m } ⊂ B for some nonnegative finite integer m and we sometimes denotes its cardinal by |x|. For further details about point processes, we refer to Daley and Vere-Jones

imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017

imsart-ejs ver. 2014/10/16 file: appIntDPP.tex date: December 1, 2017Coeurjolly and Lavancier/Intensity approximation using DPP

Acknowledgements

The authors are sincerely grateful to Adrian Baddeley and Gopalan Nair for sharing the Monte-carlo replications produced in Baddeley and Nair (2012) we used to compare the Poisson-saddlepoint approximation and the DPP approximation (Figure 3 (a),(b) and (d)). The research of J-F. Coeurjolly is supported by the Natural Sciences and Engineering Research Council of Canada.


Jerome Coste

Simone Hemm

Ulla Miguel

Philippe P Derost

François Caire

Franck Durif

Jean-Jacques Lemaire

Somatotopy in the GPi: analysis of motor side effects during intraoperative assessment in a Parkinsonian

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION MATERIAL AND METHODS

A SOMATOTOPY OF DYSTONIC MOTOR SIDE EFFECTS WAS ENCOUNTERED INSIDE EACH PALLIDAL STRUCTURE

The somatotopic organization analysis is an interesting spin-off of intra operative motor side effect analysis. The GP somatotopy related to these clinical conditions has to be confirmed. Studies have reported a somatotopy of clinical effect [dystonia:(9), dyskinesia (2)] or a spatial organization of kinesthetic cells in Parkinsonian GP (7). Techniques of localization of clinical effects or anatomical structures should influence the results and have to be refined as proposed with our patient based approach.

DISCUSSION History

A 68-year-old right-handed woman suffered from Parkinson's disease for 16 years, initially affecting the right upper limb with tremor. Before surgery, dyskinesia was encountered for the upper and the lower limbs, for the trunk with alternatively, severe blocking episodes. Moreover, an important bradykinesia was observed with chin and upper limb tremor.

Surgical technique

Stereotactic MRI: On the first day, the stereotactic Leksell G frame (Elekta, Sweden) was placed with its repositioning kit (Leksell repositioning kit, Elekta, Sweden) under local anesthesia. A stereotactic MRI was performed (Sonata 1.5 Tesla, Siemens, Germany) in three orthogonal planes (axial, coronal and sagittal). The voxel size was 0.52 × 0.62 × 2 mm 3 (field of view = 270 mm; slice thickness = 2 mm). A Turbo Spin Echo (T2-TSE) sequence was used for the axial and the sagital planes: TR = 8000 msec; TE = 10 msec; 24 images. A White matter Attenuated Inversion Recovery (WAIR) sequence was used for the coronal plane: TR = 4500 msec; TE = 13 msec; TI = 160 msec; 2 × 15 images (two interleaved series). The frame was removed.

Planning: Direct targeting of GPi was performed using a stereotactic software (Iplan, BrainLab, Germany). We identified all the segments of the Globus Pallidus (GP). GP is commonly divided in two parts (3,8) the lateral and the medial segments, respectively named external GP (GPe) and internal GP (GPi). An accessory lamina or lamina pallidi incompleta, in a basal and posterior position, partially separates the GPi into a lateral (GPil) and a medial (GPim) subdivisions (4,5,6,8). With the help of probe view reconstructions, we determined a trajectory, with a double obliquity, avoiding the main vessels, going strictly through the parenchyma and with an endpoint located at the lower boundary of the structure: we schematically placed the end point at the junction between the 3 rd and the 4 th quarter of GP according to its main antero-posterior axis (1). We planned three parallel tracts, a central one centered on the best anatomical target, a 2mm-lateral one and a 2mm-medial one on the left hemisphere, a central one, a 2mm-anterior one and a 2mm lateral one on the right hemisphere.

Implantation surgery: The following day, the stereotactic frame was repositioned. Under local anesthesia, after recording of neuronal activity, an acute stimulation (monophasic, negative current, frequency: 130Hz, pulse width: 0.06msec) was performed using unshielded tungsten microelectrodes (stimulation tip: length 1.2mm, diameter 550µm) with a microguide system (Alpha Omega, Israel). We explored the clinical effects every 1 or 2mm (current up to 5mA) along the distal 10mm on the 3 parallel tracts on each side. Controlateral dystonic movements were noted (mean current = 0.97+/-0.14mA). On the 6 tracts (45 GP checkpoints), the involved body parts (face, upper and lower limbs) were matched with the anatomic structures. Then a quadripolar electrode (DBS 3387, Medtronic, USA) was placed on the selected tract (the central one for both) with one or two contacts in the clinically most efficient area. Peroperative tele radiographic X-rays controls and postoperative non-stereotactic coronal T2-TSE MRI acquisitions were performed to control the electrode positioning.

Postoperatively, chronic DBS dramatically improved dyskinesia (electrodes implanted on the central tracks).

Surrounding of pallidal subdivisions and trajectory definition.

We found inside each structure (GPe, GPil and Post surgery clinical results.

• Initially, GPi DBS totally abolished dyskinesia.

• 3 months after the surgery ►adverse effect: dyskinesia improved by  voltage.

  lateral (GPil) and GPi medial (GPim) seemed to be characterized by a segregated body map. Controlateral dystonic movements were noted during acute intraoperative stimulation and represented as pictograms along the distal 10mm on 3 parallel tracts on each side using probe views. Central (C), anterior (A) and lateral (L) tracts are explored on the Right pallidum and central (C), medial (M) and lateral (L) trajectories on the Left one. A scale represents the 10mm course of the stimulation exploration. (figures 1b and 1c) Axial slices with surrounded structures are illustrated for the position -8mm, -4mm then 0mm on the central tracts (figures 1a and 1d).

GPim a rostro-caudal somatotopic organisation

	-face	► superior part
	-upper limb	► intermediate part
	-lower limb	► inferior part.

jcoste@chu-clermontferrand.fr References: 1.

Coubes,P. et al. 2002. Stereotact Funct Neurosurg. 4. Parent,A. 1996

. Carpenter's Human Neuroanatomy. 7. Taha,J.M. et al. 1996. Neurosurgery. 2. Kishore,A. et al. 2000. Brain. 5. Parent,M. and A.Parent. 2004. Parkinsonism Relat Disord. 8. Talairach,J. et al. 1957. Atlas d'anatomie stéréotaxique. 3. Nieuwenhuys,R. et al. 1981. The Human Central Nervous System. A synopsis and Atlas. 6. Schaltenbrand,G. and P.Bailey. 1959. Introduction to stereotaxis with an atlas of the human brain. 9. Vayssiere,N. et al. 2004. J Neurosurg.


Florian Huet
email: florian.huet@univ-smb.fr

Émile Roux

Aurélien Carré

Small nomadic counter-rotating wind turbine : Design and experimentation

Keywords: Small-scale wind turbine, counter-rotating propeller, wind energy, nomadic power

Many studies show the efficiency of large-scale wind turbines in wind energy harvesting for collective applications. However, few examples relate to small one stage wind turbines (OSWT) for mobile and individual applications. This work details the modeling and the prototype realization of a 18 cm diameter wind turbine using counter-rotating propeller technology (CRWT). This device has been tested in a wind tunnel with wind velocities in the range of 3.5 to 8 m.s -1 . First experimental results with a wind of 4.2 m.s -1 (which is the theoretical optimum wind velocity of the device) show that the wind turbine can provide a power of 0.6 W currently dissipated in a load resistance of 1.5 Ω. The power goes up to 2.8 W for a wind of 8 m.s -1 .

I. INTRODUCTION

The green sources of energy are a major concern and become even a necessity while we wish to deploy an activity in nomadic or isolated natural environment. The energy demand goes from some watts to recharge batteries of low capacity (smartphone, GPS, lamp…) to some hundred watts to feed isolated structures as a mountain refuge for example. The technologies usually used for these applications exploit the solar or wind energy. The power density supplied by the outdoor solar energy can be very important, 15 000 µW.cm -2 , but it decreases strongly in 150 µW.cm -2 when weather conditions are degraded [START_REF] Roundy | A study of low level vibrations as a power source for wireless sensor nodes[END_REF]. On the other hand, the power density supplied by a laminar airflow is stable is estimated at 380 µW.cm -3 for a wind velocity of 5 m.s -1 . The ubiquity of the natural wind in the relief or the sea regions favors the exploitation of this energy source. This natural wind velocity follows a Weibull distribution with average value is 3.94 m.s -1 and the maximal is about 5 m.s -1 . The source of airflow can also come from human motions and thus capture relative flows. In these two cases, the harvesting energy depends strongly on the velocity of the source and less on the apparent surface.

Many studies show the efficiency of large-scale wind turbines in wind energy harvesting for collective applications. However, few examples [START_REF] Kishore | Design and experimental verification of a high efficiency small wind energy portable turbine (SWEPT)[END_REF] relate to small wind turbines for mobile and individual applications. This paper suggests a mechanical concept of counter-rotating propeller to a small-scale wind turbine (< 20 cm). The counter-rotating propeller solution is well known in the aviation to increase the thrust of aircrafts for the same initial mechanical power level. The efficiency of this propeller is improved by 30 % compared with a conventional propeller [START_REF] Davenport | Analysis of counter-rotating propeller performance[END_REF]. The integration of this structure type of large- scale wind turbine was already realized to validate this theory [START_REF] Debleser | Wind turbine with counter rotating rotors[END_REF][START_REF] Kumar | Computational analysis of 30 kW contra rotor wind turbine[END_REF][START_REF] Hoang | An evaluation of the performance of 10kw counter-rotating wind turbine using CFD simulation[END_REF][START_REF] Mitulet | Wind Tunnel Testing for a New Experimental Model of Counter-Rotating Wind Turbine[END_REF]. However, the use of this technology on small-scale devices is a real innovation and allows envisaging an efficiency improvement of the same factor.

II. COUNTER-ROTATING WIND TURBINE

A. Counter-rotating wind turbine concept

The CRWT concept exploits a general technology adapted to propulsion propellers in aviation. In this paper, this structure is composed of two different and optimized propeller stages which rotate in the opposite direction (Figure 1a). The stages directly engrave the two independent generators. The incoming air flow is partially captured by the first propeller stage and is converted into mechanical energy by the shaft rotation. The outgoing airflow is deflected, creating a flow composed of an axial and tangential velocity. The second propeller stage is adapted to harvest this new form of wind field, but it turns in the opposite direction. The optimization of the geometry of the propeller blades maximizes the caption of wind kinetic energy. The maximum power is theoretically restricted by the Betz limit and it can be formulated from the following relation:

(1)

B. Propellers modelisation

The propellers geometries are optimized with Heliciel software. The latter proposes a reverse design of the blades in order to optimize the mechanical power according to a defined wind field (axial and tangential velocities). The characteristics of the upstream and downstream propeller stages are listed in Table 1. After the first stage incidence angle optimization, the airflow data output is injected in the second stage input data. The axial velocity component is greatly reduced compared to the source. The global geometrical parameters listed in the Table 1 are determined using the design of experiments method (DOE).

C. Electromecanical conversion

For electromechanical conversion, two independent permanent magnet synchronous generators with an optimal rotation velocity of 7200 rpm are used. The frequency of the electrical signals informs directly about the rotation velocity of the propellers. The generator is modeled (Figure 2) by three coils in series with three resistors. To harvest the electrical energy produced, a passive extraction circuit composed of a diode bridge suitable for synchronous generators is envisaged. But actually, the electrical voltages at the coil terminals are too low to pass the threshold diode voltage. The power is measured on an optimum load resistance of 1.5 Ω in parallel with each coil.

III. EXPERIMENTAL RESULTS

The tests are carried out in a wind tunnel with an air stream of 23×23 cm 2 . An OSWT and a CRWT (Figure 1b) prototype were subjected to a wind source ranging from 3.5 m.s -1 to 8 m.s -1 . The Figure 3 shows the produced electrical power by the OSWT and the optimized CRWT and the maximum recoverable power P max (Equation 1) as a function of the wind velocity. The classical wind turbine starts at 4.6 m.s -1 and it produce 0.01 W. The counter-rotating starts around 4.2 m.s -1 with small difference between the two stages. This device produces 0.6 W which corresponds to the maximum efficiency of 85 %. At 8 m.s -1 , the OSWT produces 1.3 W and the counter-rotating produces 2.8 W.

IV. CONCLUSIONS AND PROSPECTS

A prototype of a counter-rotating wind turbine was realized and tested. The latter is above all a means of observing and measuring the effects of the airflow on the two stages independently. The first experimental test show promising results with an electrical power of 0.6 W for a wind velocity of 4.2 m.s -1 and a maximum power of 2.8 W at the maximum test velocity. The optimum efficiency is 85 %. A semi-analytical model integrating the optimized propellers theoretical data is under development, in order to design a new small dimension counter-rotating wind turbine. These first results also show the need to integrate an intermediate transmission composed of a mechanical amplifier and a mechanical differential to couple both stages. This step will be included in the next prototype in order to exploit optimally a single generator while dissociating the propeller stages.

FIGURE 1 .

 1 FIGURE 1. COUNTER-ROTATING WIND TURBINE CONCEPT AND PROTOTYPE

FIGURE 3 .

 3 FIGURE 3. EXPERIMENTAL POWER RESULTS

TABLE 1 .

 1 PROPELLERS CHARACTERISTICS

		ROPELLER STAGE
	PARAMETERS	UPSTREAM	DOWNSTREAM
	BLADE NUMBER	6	6
	PROFIL	NACA1412	NACA6412
	BLADE FOOT / HEAD DIAMETER	40 mm / 180 mm
	BLADE FOOT / HEAD WIDTH	17 mm / 7 mm	17 mm / 17 mm
	THEORETICAL MAXIMUM POWER	0.8 W	2.2 W
	THEORETICAL ROTATION VELOCITY	2800 rpm	6000 rpm
	FIGURE 2.

WIND TURBINE ELECTRICAL SCHEMATISATION

ACKNOWLEDGMENT

The authors would like to thank L. Bouderaux, A. Brenin, T. Ringeisen and L. Roques, students of Polytech Annecy, for their contributions in INoWind project and the MPH department of Annecy IUT for the wind tunnel provision.


Florian Huet

Emile Roux

Aurélien Carré

Aurélien Carré Small

Small nomadic counter-rotating wind turbine : design and experimentation

come

1 Contact : Florian.Huet@univ-smb.fr

  Small nomadic counter-rotating wind turbine : Design and experimentation Florian Huet 1 , Emile Roux and Aurélien Carré Univ. Savoie Mont Blanc, SYMME, F74000, Annecy, a new high efficiency small scale wind turbine generator providing a power of a few Watts power Application : Recharge low capacity batteries (Smartphone, GPS, lamp…) authors would like to thank L. Bouderaux, A. Brenin, T. Ringeisen and L. Roques, students of Polytech Annecy, for their contributions in INoWind project and the MPH department of Annecy IUT for the wind tunnel provision.

First

  One Stage Wind Turbine (OSWT) :𝑃 𝑂𝑆𝑊𝑇 = 0 W at 4,2 m.s -1 , 1.3 W at 8 m.s -1Counter-rotating Wind Turbine (CRWT) :𝑃 𝐶𝑅𝑊𝑇 = 0.6 W at 4,2 m.s -1 , 2.8 W at 8 m.s -1 Maximum efficiency 85 % at 4,2 m.s -1 Speed : Integrate a venturi into the prototype to increase wind speed and concentrate the air flow Improve The Transmission : Use a single generator while separating both stages thanks to a mechanical differential