Mariia Vladimirova
email: mariia.vladimirova@inria.fr

Stéphane Girard

Hien Nguyen

Julyan Arbel

Sub-Weibull distributions: generalizing sub-Gaussian and sub-Exponential properties to heavier-tailed distributions

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction and definition

Sub-Gaussian distributions, respectively sub-Exponential, are characterized by their tails being upper bounded by Gaussian, respectively Exponential, tails. More precisely, we say that a random variable X is sub-Gaussian, resp. sub-Exponential, if there exist positive constants a and b such that P(|X| ≥ x) ≤ a exp(-bx 2 ), resp. P(|X| ≥ x) ≤ a exp(-bx), for all x > 0.

(1)

These properties have been intensely studied in the recent years due to their relationship with various fields of probability and statistics, including concentration, transportation and PAC-Bayes inequalities [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF][START_REF] Raginsky | Concentration of measure inequalities in information theory, communications, and coding[END_REF][START_REF] Van Handel | Probability in high dimension[END_REF][START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF], the missing mass problem [START_REF] Ben-Hamou | Concentration inequalities in the infinite urn scheme for occupancy counts and the missing mass, with applications[END_REF], bandit problems [START_REF] Bubeck | Regret Analysis of Stochastic and Nonstochastic Multi-Armed Bandit Problems[END_REF] and singular values of random matrices [START_REF] Rudelson | Non-asymptotic theory of random matrices: extreme singular values[END_REF].

It is tempting to generalise (1) by considering the class of distributions satisfying P(|X| ≥ x) ≤ a exp -bx 1/θ , for all x > 0, for some θ, a, b > 0.

(2) This is the goal of the present note. Since a Weibull random variable X on R + is defined by a survival function, for x > 0, F (x) = P(X ≥ x) = exp -bx 1/θ , for some b > 0, θ > 0,

we term a distribution satisfying (2) a sub-Weibull distribution (see [START_REF] Rinne | The Weibull Distribution: a Handbook[END_REF], for a detailed account on the Weibull distribution). An illustration of symmetric sub-Weibull distributions is represented in Figure 1 for different values of the tail parameter θ. Since only the tail is relevant for the illustration, the sub-Weibull survival functions S θ depicted here are obtained from the survival functions S W θ of Weibull random variables with shape parameter 1/θ, scale parameter 1. More specifically, only the part of S W θ to the right of the 95-th quantile is considered, and it is symmetrized around zero:

S θ (x) = 10e -(x+log(20) θ ) 1/θ if x ≥ 0, and S θ (x) = 1 -10e -(-x+log(20) θ ) 1/θ if x < 0.
Proof. In the proof we use the notation between two positive sequences (a k ) k and (b k ) k , writing a k b k , if there exists a constant C > 0 such that for all integer k, a k ≤ Cb k . 1 ⇒ 2. Assume Property 1 holds. We have, for any k ≥ 1,

E |X| k (a) = ∞ 0 P |X| k > x dx = ∞ 0 P |X| > x 1/k dx (b) ≤ ∞ 0 2 exp -(x/K 1 ) 1/(kθ) dx = 2K 1 kθ ∞ 0
e -u u kθ-1 du = 2K 1 kθ Γ (kθ) = 2K 1 Γ (kθ + 1)

(c)
2K 1 (kθ + 1) kθ+1 = 2K 1 (kθ + 1) (kθ + 1) kθ , where (a) is by the so-called integral identity for |X| k , (b) is by Property 1, and (c) comes from Stirling's formula, yielding Γ(u) u u . Taking the k-th root of the expression above yields

X k (2K 1 ) 1/k (kθ + 1) 1/k (kθ + 1) θ k θ ,
which is to Property 2. 2 ⇒ 3. Assume Property 2 holds. Recalling the Taylor series expansion of the exponential function, we obtain

E exp λ 1/θ |X| 1/θ = E 1 + ∞ k=1 (λ 1/θ |X| 1/θ ) k k! = 1 + ∞ k=1 λ k/θ E[|X| k/θ ] k! . Property 2 guarantees that E[|X| k/θ ] ≤ K k/θ
2 (k/θ) k for some K 2 and for any k ≥ θ. There exist only a finite number of integers k that are less than θ, thus there exists a constant K2 such that

E[|X| k/θ ] ≤ Kk/θ 2 (k/θ) k for integers k ≥ θ. Defining K 2 = max{K 2 , K2 } implies E[|X| k/θ ] ≤ K k/θ 2
(k/θ) k for all integers k ≥ 1. By Stirling's approximation we have k! (k/e) k . Substituting these two bounds, we get

E exp λ 1/θ |X| 1/θ 1 + ∞ k=1 λ k/θ K k/θ 2 (k/θ) k (k/e) k = ∞ k=0 K k/θ 2 (eλ 1/θ /θ) k = 1 1 -e(K 2 λ) 1/θ /θ
provided that e(K 2 λ)1/θ /θ < 1, in which case the geometric series above converges. To bound this quantity further, we can use the inequality 1 1-x ≤ e 2x , which is valid for x ∈ [0, 1/2]. It follows that

E exp λ 1/θ X 1/θ ≤ exp 2e(K 2 λ) 1/θ /θ for all λ satisfying 0 < λ ≤ 1 K 2 θ 2e
θ and some positive K 2 defined above. This yields Property 3 with

K 3 = K 2 (2e/θ) θ .
3 ⇒ 4. Assume Property 3 holds. Take λ = 1/K 4 , where K 4 ≥ K 3 /(ln 2) θ . This yields Property 4.

4 ⇒ 1. Assume Property 4 holds. Then, by Markov's inequality and Property 3, we obtain

P |X| > x = P e (|X|/K4) 1/θ > e (x/K4) 1/θ ≤ E[e (|X|/K4) 1/θ ] e (x/K4) 1/θ ≤ 2e -(x/K4) 1/θ .
This proves Property 1 with

K 1 = K 4 .
Informally, the tails of a subW (θ) distribution are dominated by (i.e. decay at least as fast as) the tails of a Weibull variable with shape parameter 1 equal to 1/θ. Sub-Gaussian and sub-Exponential variables, which are commonly used, are special cases of sub-Weibull random variables with tail parameter θ = 1/2 and θ = 1, respectively, see Table 1. Sub-Gaussian distributions are Table 1: Sub-Gaussian, sub-Exponential and sub-Weibull distributions comparison in terms of tail P (|X| ≥ x) and moment condition, with K 1 and K 2 some positive constants. The first two are a special case of the last with θ = 1/2 and θ = 1, respectively.

Distribution

Tails Moments

Sub-Gaussian

P(|X| ≥ x) ≤ 2e -(x/K1) 2 X k ≤ K 2 √ k Sub-Exponential P(|X| ≥ x) ≤ 2e -x/K1 X k ≤ K 2 k Sub-Weibull P(|X| ≥ x) ≤ 2e -(x/K1) 1/θ X k ≤ K 2 k θ
sub-Exponential as well. Such inclusion properties are generalized to sub-Weibull distributions in the following proposition.

Proposition 2.1 (Inclusion). Let θ 1 and θ 2 such that 0 < θ 1 ≤ θ 2 be two sub-Weibull tail parameters. The following inclusion holds subW (θ 1 ) ⊂ subW (θ 2 ).

Proof. For X ∼ subW (θ 1 ), there exists some constant K 2 > 0 such that for all k > 0,

X k ≤ K 2 k θ1 . Since k θ1 ≤ k θ2 for all k ≥ 1, this yields X k ≤ K 2 k θ2 , which by definition implies X ∼ subW (θ 2 ).
Let a random variable X follow a sub-Weibull distribution with tail parameter θ. Due to the property of inclusion from Proposition 2.1, the sub-Weibull definition provides an upper bound for the tail. In order to address the question of a lower bound on the tail parameter of some sub-Weibull distribution, we rely on an optimal tail parameter for that distribution through the use of the moment Property 2 of Theorem 2.1. To this aim, we introduce the notation of asymptotic equivalence between two positive sequences (a k ) k and (b k ) k as by a k b k , defined by the existence of constants

D ≥ d > 0 such that d ≤ a k b k ≤ D, for all k ∈ N. (4)
We can now introduce the notion of optimal sub-Weibull tail coefficient, and provide a momentbased condition for optimality to hold.

Proposition 2.2 (Optimal sub-Weibull tail coefficient and moment condition). Let θ > 0 and let X be a random variable satisfying the following asymptotic equivalence on moments

X k k θ .
Then X is sub-Weibull distributed with optimal tail parameter θ, in the sense that for any θ < θ, X is not subW (θ ).

Proof. By the upper bound of the asymptotic equivalence assumption on moments, X satisfies Property 2 of Theorem 2.1, so X ∼ subW (θ). Let θ < θ. By the lower bound of the asymptotic equivalence assumption on moments, there does not exist any constant K 2 such that X k ≤ K 2 k θ for any k ∈ R, so X is not sub-Weibull with tail parameter θ . This concludes the proof.

When sub-Gaussian and sub-Exponential properties are studied, it is often assumed that the distribution is centered. However, note that these properties, as well as the sub-Weibull property, are tail properties, and as such, they are not affected by translation. Non-centered random variables can always be centered by subtracting the mean without affecting their tail property. In particular, variable centering does not change the optimal tail parameter of a sub-Weibull distribution. Scaling procedure does not influence the tail neither: multiplication by some fixed number will change the coefficient K 1 in the tail Property 1 of Theorem 2.1, but not the tail parameter.

It is known that the product of sub-Gaussian random variables is sub-exponential (cf. Vershynin 2018, Lemma 2.7.7). By the same virtue, we have the fact that the class of sub-Weibull random variables is closed under multiplication and addition.

Proposition 2.3 (Closure under multiplication and addition). Let X and Y be sub-Weibull random variables with tail parameters θ 1 and θ 2 , respectively. Then, XY and X + Y are sub-Weibull with respective tail parameters θ 1 + θ 2 and max(θ 1 , θ 2 ).

Proof. Let us focus on the closure under addition property. Property 2 in Theorem 2.1 and the triangular inequality entail that there exist K 2,1 > 0 and K 2,2 > 0 such that, for all k ≥ 1, θ1,θ2) .

X + Y k ≤ X k + Y k ≤ K 2,1 k θ1 + K 2,2 k θ2 ≤ (K 2,1 + K 2,2 )k max(
The conclusion follows from Theorem 2.1. Furthermore, we may establish that the class of sub-Weibull random variables is larger than the class of bounded random variables but is smaller than the class of random variables with finite p-th moment, for every p ∈ [0, ∞). This result is comparable to Remark 2.7.14 of [START_REF] Vershynin | High-Dimensional Probability: An Introduction with Applications in Data Science[END_REF], which establishes the same relationship regarding the class of sub-Gaussian random variables.

Concentration properties

In this section, we focus on theoretical results about the sum of sub-Weibull random variables, including concentration properties.

Proposition 3.1 (Sum of sub-Weibull random variables). Let X 1 , . . . , X n be sub-Weibull random variables with tail parameter θ. Then the sum n i=1 X i is sub-Weibull with tail parameter θ. Proof. Using Property 2 in Theorem 2.1 combined with Minkowski's inequality yields (for any k ≥ 1):

n i=1 X i k ≤ n i=1 X i k ≤ n i=1 K 2,i k θ = K 2,Σ k θ , (5)
where K 2,Σ = n i=1 K 2,i . It implies that the sum X 1 + • • • + X n satisfies the sub-Weibull property with tail parameter at most θ. In addition, if X 1 , . . . , X n are from the same distribution, i.e. X i k ≤ K 2 k θ for all i ∈ {1, . . . , n} with some constant K 2 > 0, then K 2,Σ = nK 2 . Using Proposition 2.3, and the proof of Property 2 in Theorem 2.1, we have the following Hoeffding-type concentration inequality regarding the sum n i=1 X i .

Corollary 3.1 (Concentration of the sum). Let that X 1 , . . . , X n be identically distributed sub-Weibull random variables with tail parameter θ. Then, for all x ≥ nK θ , we have

P n i=1 X i ≥ x ≤ exp - x nK θ 1/θ
, for some constant K θ dependent on θ.

Proof. For all real k ≥ 1, the bound (5) in the proof of Proposition 3.1 yields

E   n i=1 X i k   ≤ K 2,Σ k θ k = nK 2 k θ k ,
and thus, by Markov's inequality,

P n i=1 X i ≥ t = P   n i=1 X i k ≥ t k   ≤ E | n i=1 X i | k t k ≤ (nK 2 ) k k θ k t k .
Choose t so that exp (-k) = (nK 2 ) k k θ k /t k . Then, we have

P n i=1 X i ≥ enK 2 k θ ≤ exp (-k) . ( 6
) Let K θ = eK 2 and notice that nK θ k θ = x implies that k = x nK θ 1/θ
and thus, for x ≥ nK θ (since k ≥ 1):

P n i=1 X i ≥ x ≤ exp - x nK θ 1/θ
, which is the expected result.

Let us remark that, alternatively, letting exp (-k) = α in (6) implies, for 0 < α < 1/e, the confidence statement

P n i=1 X i ≤ nK θ log 1 α θ ≥ 1 -α.
Under the same conditions as Proposition 3.1, under the restriction that θ ≤ 1, [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF] proposed the inequality

P n i=1 X i ≥ x ≤ K θ exp - 1 K θ min x 2 n , x 1/θ n (1-θ)/θ , ( 7
)
where K θ is again a constant that only depends on θ. We observe that (1θ) /θ < 1/θ and thus the bound of (7) converges faster to zero than that of Corollary 3.1 in all cases where they are comparable. A version of this result also appears in Lemma 3.5 of [START_REF] Adamczak | Restricted isometry property of matrices with independent columns and neighborly polytopes by random sampling[END_REF]. Further results regarding the concentration of the sum of sub-Weibull random variables are available in [START_REF] Kuchibhotla | Moving beyond sub-Gaussianity in highdimensional statistics: Applications in covariance estimation and linear regression[END_REF] and [START_REF] Hao | Bootstrapping upper confidence bound[END_REF]. See also [START_REF] Bakhshizadeh | Sharp concentration results for heavy-tailed distributions[END_REF] for further concentration results regarding heavy-tailed distributions.

Application to Bayesian neural networks

This section gives an example of sub-Weibull variables that arise in deep learning, more specifically in the context of Bayesian deep neural networks, as originally described in [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF]. According to Theorem 2.1 of [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF], they correspond respectively to subW ( 1 /2), subW (1), subW ( 3 /2), subW (5) and subW (50).

Distribution at the level of units

We first describe so-called fully connected neural networks. Neural networks are hierarchical models made of layers: an input, several hidden layers and an output. Each layer following the input layer consists of units which are linear combinations of previous layer units transformed by a nonlinear function, often referred to as the nonlinearity or activation function denoted by φ : R → R. Given an input x ∈ R N (for instance an image made of N pixels) the -th hidden layer consists of two vectors whose size is called the width of layer, denoted by H . The vector of units before application of the non-linearity is called pre-nonlinearity, and is denoted by g ( ) = g ( ) (x), while the vector obtained after element-wise application of φ is called post-nonlinearity and is denoted by h ( ) = h ( ) (x). More specifically, these vectors are defined as

g ( ) (x) = W ( ) h ( -1) (x), h ( ) (x) = φ(g ( ) (x)), (8)
where W ( ) is a weight matrix of dimension H × H -1 including a bias vector, with the convention that H 0 = N , the input dimension.

Neural networks perform state-of-the-art results in many areas. Researchers aim at better understanding the driving mechanisms behind their effectiveness. In particular, the study of the neural networks distributional properties through Bayesian analysis, where weights are assumed to follow a prior distribution, has attracted a lot of attention in the recent years. The main focus of research in the field shows that a Bayesian deep neural network converges in distribution to a Gaussian process when the width of all the layers goes to infinity. See for instance [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF]; [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF] for the main proofs. Further research such as [START_REF] Hayou | On the impact of the activation function on deep neural networks training[END_REF] builds upon the limiting Gaussian process property of neural networks, as well as on the notion of edge of chaos developed by [START_REF] Schoenholz | Deep information propagation[END_REF], in order to devise novel architecture rules for neural networks.

In contrast, Theorem 2.1 of [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF] shows that the non asymptotic (i.e. for finite width neural networks) prior distribution of units from the -th layer (both before and after activation, g ( ) and h ( ) ) induced by a standard Gaussian prior on the weights W ( ) , is sub-Weibull with tail parameter /2. Therefore, the deeper the layer is, the heavier-tailed is the units distribution. This result puts into perspective the infinite width Gaussian process property which might be far from holding for real world neural networks.

In order to illustrate the sub-Weibull property of units prior distributions, we performed the following experiment with a deep neural network of 100 layers. We considered a Bayesian neural network with independent standard Gaussian priors on the weights, where layers are composed of H = 1000 -10( -1) hidden units, = 1, . . . , 100, with the ReLU nonlinearity φ defined by φ(x) = max(0, x) (see ( 8)). The input vector x contains 10 4 numbers sampled from independent standard Gaussian. In order to evaluate the units prior distributions, we used a Monte Carlo approximation, where the input vector x was kept fixed, and the weights were sampled from independent standard Gaussian n = 10 5 times. Figure 2 illustrates the survival function of prenonlinearity hidden units g ( ) for layers = 1, 2, 3, 10 and 100. This indicates that the prior tails of the units get heavier when they originate from deeper layers, in accordance with the main result of [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF].

Tail parameter estimation

We conclude this section by suggesting a statistical estimation procedure for the tail parameter θ. We adopt an approach that relies on Weibull random variables with tail parameter θ > 0 and scale parameter λ > 0, which are subW (θ). The cumulative distribution function is F (y) = 1-e -(y/λ) 1/θ for all y ≥ 0, with corresponding quantile function q(t) = λ (-log(1t))

θ . Taking the logarithm, we obtain the following expression

log q(t) = θ log log 1 1 -t + log λ,
showing an affine relationship involving the log quantiles and θ as slope parameter. This suggests a simple estimation procedure for θ by linear regression. More specifically, assume we have n iid observations Y 1 , . . . , Y n , and denote the order statistics by Y i,n . Consider a fraction k ≤ n of the largest observations, that are the order statistics Y n-i+1,n for 1 ≤ i ≤ k. These are approximating the quantiles of order n-i n = 1 -i n , respectively. As a consequence, we consider the slope coefficient in the regression of log Y n-i+1,n against log log(n/i), for 1 ≤ i ≤ k, as an estimate of θ. This estimator was first introduced by [START_REF] Gardes | Estimation of the Weibull tail-coefficient with linear combination of upper order statistics[END_REF] in the context of Weibull-tail distributions. We visualize the alignment in a quantile-quantile plot in log-log scale in Figure 3.

We considered the following experiment in order to illustrate the estimation procedure. We have built neural networks of 10 hidden layers, with H = H hidden units on each layers, and made H vary in {1, 3, 10}. We used a fixed input x of size 10 4 , which can be thought of as an image of dimension 100 × 100. This input was sampled once for all with standard Gaussian entries. In order to obtain samples from the prior distribution of the neural network units, we have sampled the weights from independent centered Gaussians with variance set to 2, from which units were obtained by forward evaluation with the ReLU non-linearity. This process was iterated n = 10 5 times and the number of largest observations selected is k = 10 3 . We show in Figure 3 the θ estimates obtained by using pre-nonlinearities g ( ) of layers in {1, 3, 10} as data, for H varying in {1, 3, 10}. We can see that the theoretical result from [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF] that states that θ = /2 is well in line with the estimates obtained with networks of width H = 1. When the width increases, the estimates for θ tend to decrease, narrowing the gap to the lower bound of 1/2. These results illustrate the large width results by [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF]; [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF], which state that in the large width regime, the unit distribution tends to a Gaussian, thus with a limiting θ parameter of 1/2. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.6

1.8 2.0 2.2 2.4 1.0 1.1 1.2 1.3 1.4 l = 1, H = 1 θ ^ = 0.58 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.6 1.8 2.0 2.2 2.4

1.5 2.0 2.5 3.0 l = 3, H = 1 θ ^ = 1.71 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.6

1.8 2.0 2.2 2.4 1 2 3 4 5 l = 10, H = 1 θ ^ = 4.37 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 1.2 1.4 1.6 1.8 2.0 2.2 l = 1, H = 3 θ ^ = 0.66 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 4.0 4.5 5.0 l = 3, H = 3 θ ^ = 1.18 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 12.0 13.0 14.0 l = 10, H = 3 θ ^ = 2.16 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 1.2 1.4 1.6 1.8 2.0 l = 1, H = 10 θ ^ = 0.67 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 5.0 5.4 5.8 l = 3, H = 10 θ ^ = 0.89 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 Figure 3: Illustration of the estimation procedure for the tail parameter θ for the prior distribution of Bayesian neural networks units. Layers of level 1, 3 and 10 on the left, middle and right panels, respectively. Width H (for all considered layers ) equal to 1, 3 and 10 on the top, middle and bottom rows, respectively. Note that the theoretical tail parameter value stated by [START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF] is θ = /2, while the infinite width (H → ∞) limiting value is θ = 1/2.

Figure 1 :

 1 Figure 1: Illustration of sub-Weibull survival curves on R with varying tail parameters θ.

Figure 2 :

 2 Figure2: Illustration of layers = 1, 2, 3, 10 and 100 hidden units marginal prior distributions tails. According to Theorem 2.1 of[START_REF] Vladimirova | Understanding priors in Bayesian neural networks at the unit level[END_REF], they correspond respectively to subW ( 1 /2), subW (1), subW ( 3 /2), subW (5) and subW (50).

Weibull distributions are commonly parameterized by a shape parameter κ. Here we use instead θ = 1/κ for the convenience that the larger the tail parameter θ, the heavier the tails of the sub-Weibull distribution.

are funded by Grenoble Alpes Data Institute, supported by the French National Research Agency under the "Investissements d'avenir" program (ANR-15-IDEX-02). H. Nguyen is funded by the Australian Research Council grants: DE170101134 and DP180101192.


T C Phi

A Muzy

P Reynaud-Bouret

Event-scheduling algorithms with Kalikow decomposition for simulating potentially infinite neuronal networks

Keywords: MSC 2010 subject classifications:, Kalikow decomposition, Discrete event simulation, Point process, Infinite neuronal networks

, which always needs to simulate the whole network to access the behaviour of one particular neuron of the network. On the other hand, for discrete time models, theoretical algorithms based on Kalikow decomposition can pick at random influencing neurons and perform a perfect simulation (meaning without approximations) of the behaviour of one given neuron embedded in an infinite network, at every time step. These algorithms are currently not computationally tractable in continuous time. To solve this problem, an event-scheduling algorithm with Kalikow decomposition is proposed here for the sequential simulation of point processes neuronal models satisfying this decomposition. This new algorithm is applied to infinite neuronal networks whose finite time simulation is a prerequisite to realistic brain modeling.

Introduction

Point processes in time are stochastic objects that model efficiently event occurrences with a huge variety of applications: time of deaths, earthquake occurrences, gene positions on DNA strand, etc. [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF][START_REF] Vere-Jones | Some examples of statistical estimation applied to earthquake data[END_REF][START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF]).

Most of the time, point processes are multivariate [START_REF] Didelez | Graphical models of markes point processes based on local independence[END_REF] in the sense that either several processes are considered at the same time, or in the sense that one process regroups together all the events of the different processes and marks them by their type. A typical example consists in considering either two processes, one counting the wedding events of a given person and one counting the children birth dates of the same person or only one marked process which regroups all the possible dates of birth or weddings independently and adds one mark per point, here wedding or birth.

Consider now a network of neurons each of them emitting action potentials (spikes). These spike trains can be modeled by a multivariate point process with a potentially infinite number of marks, each mark representing one given neuron. The main difference between classical models of multivariate point processes and the ones considered in particular for neuronal networks is the size of the network. A human brain consists in about 10 11 neurons whereas a cockroach contains already about 10 6 neurons. Therefore the simulation of the whole network is either impossible or a very difficult and computationally intensive task for which particular tricks depending on the shape of the network or the point processes have to be used [START_REF] Peters | Rejection-free MonteCarlo sampling for general potentials[END_REF][START_REF] Dassios | Exact simulation of Hawkes process with exponentially decaying intensity[END_REF][START_REF] Mascart | Centralized and distributed simulations of point processes using local independence graphs: A computational complexity analysis[END_REF].

Another point of view, which is the one considered here, is to simulate, not the whole network, but the events of one particular node or neuron, embedded in and interacting with the whole network. In this sense, one might consider an infinite network. This is the mathematical point of view considered in a series of papers [START_REF] Galves | Infinite Systems of Interacting Chains with Memory of Variable Length A Stochastic Model for Biological Neural Nets[END_REF][START_REF] Galves | Modeling networks of spiking neurons as interacting processes with memory of variable length[END_REF][START_REF] Ost | Sparse space-time models: concentration inequalities and Lasso[END_REF] and based on Kalikow decomposition [START_REF] Kalikow | Random markov processes and uniform martingales[END_REF] coupled with perfect simulation theoretical algorithms [START_REF] Comets | Processes with long memory: Regenerative construction and perfect simulation[END_REF][START_REF] Fernandéz | Coupling, renewal and perfect simulation of chains of infinite order[END_REF]. However these works are suitable in discrete time and only provide a way to decide at each time step if the neuron is spiking or not. They cannot operate in continuous time, i.e. they cannot directly predict the next event (or spike). Up to our knowledge, there exists only one attempt of using such decomposition in continuous time [START_REF] Hodara | Hawkes Processes with variable length memory and an infinite number of components[END_REF], but the corresponding simulation algorithm is purely theoretical in the sense that the corresponding conditional Kalikow decomposition should exist given the whole infinite realization of a multivariate Poisson process, with an infinite number of marks, quantity which is impossible to simulate in practice.

The aim of the present work is to present an algorithm which

• can operate in continuous time in the sense that it can predict the occurrence of the next event. In this sense, it is an event-scheduling algorithm; • can simulate the behavior of one particular neuron embedded in a potentially infinite network without having to simulate the whole network; • is based on an unconditional Kalikow decomposition and in this sense, can only work for point processes with this decomposition.

In Section 2, we specify the links between event-scheduling algorithms and the classical point process theory. In Section 3, we give the Kalikow decomposition. In Section 4, we present the backwardforward perfect simulation algorithm and prove why it almost surely ends under certain conditions. In Section 5, we provide simulation results and a conclusion is given in Section 6.

Event-scheduling simulation of point processes

On the one hand, simulation algorithms of multivariate point processes [START_REF] Ogata | On Lewis Simulation Method for Point Processes[END_REF] are quite well known in the statistical community but as far as we know quite confidential in the simulation (computer scientist) community. On the other hand, event-scheduling simulation first appeared in the mid-1960s [START_REF] Tocher | PLUS/GPS III Specification[END_REF] and was formalized as discrete event systems in the mid-1970s [START_REF] Zeigler | Theory of Modelling and Simulation[END_REF] to interpret very general simulation algorithms scheduling "next events". A basic event-scheduling algorithm "jumps" from one event occurring at a time stamp t ∈ R + 0 to a next event occurring at a next time stamp t ′ ∈ R + 0 , with t ′ ≥ t. In a discrete event system, the state of the system is considered as changing at times t, t ′ and conversely unchanging in between [START_REF] Zeigler | Theory of Modeling and Simulation: Discrete Event & Iterative System Computational Foundations[END_REF]. In [START_REF] Mascart | Centralized and distributed simulations of point processes using local independence graphs: A computational complexity analysis[END_REF], we have written the main equivalence between the point processes simulation algorithms and the discrete event simulation set-up, which led us to a significant improvement in terms of computational time when huge but finite networks are into play. Usual event-scheduling simulation algorithms have been developed considering independently the components (nodes) of a system. Our approach considers new algorithms for activity tracking simulation [START_REF] Muzy | Exploiting activity for the modeling and simulation of dynamics and learning processes in hierarchical (neurocognitive) systems[END_REF]. The event activity is tracked from active nodes to children (influencees).

Here we just recall the main ingredients that are useful for the sequel.

To define point processes, we need a filtration or history (F t ) t≥0 . Most of the time, and this will be the case here, this filtration (or history) is restricted to the internal history of the multivariate process (F int t ) t≥0 , which means that at time t-, i.e. just before time t, we only have access to the events that have already occurred in the past strictly before time t, in the whole network. The conditional intensity, φ i (t|F int t-), of the point process, representing neuron i gives the instantaneous firing rate, that is the frequency of spikes, given the past contained in F int t-. Let us just mention two very famous examples.

If φ i (t|F int t-) is a deterministic constant, say M , then the spikes of neuron i form a homogeneous Poisson process with intensity M . The occurrence of spikes are completely independent from what occurs elsewhere in the network and from the previous occurrences of spikes of neuron i.

If we denote by I the set of neurons, we can also envision the following form for the conditional intensity:

φ i (t|F int t-) = ν i + j∈I w ij (Nb j [t-A,t) ∧ M ). (2.1)
This is a particular case of generalized Hawkes processes [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF]. More precisely ν i is the spontaneous rate (assumed to be less than the deterministic upper bound M > 1) of neuron i. Then every neuron in the network can excite neuron i: more precisely, one counts the number of spikes that have been produced by neuron j just before t, in a window of length A, this is Nb j [t-A,t) ; we clip it by the upper bound M and modulate its contribution to the intensity by the positive synaptic weight between neuron i and neuron j, w ij . For instance, if there is only one spike in the whole network just before time t, and if this happens on neuron j, then the intensity for neuron i becomes ν i + w ij . The sum over all neurons j mimics the synaptic integration that takes place at neuron i. As a renormalization constraint, we assume that sup i∈I j∈I w ij < 1. This ensures in particular that such a process has always a conditional intensity bounded by 2M .

Hence, starting for instance at time t, and given the whole past, one can compute the next event in the network by computing for each node of the network the next event in absence of other spike apparition. To do so, remark that in absence of other spike apparition, the quantity φ i (s|F int s-) for s > t becomes for instance in the previous example

φ abs i (s, t) = ν i + j∈I w ij (Nb j [s-A,t) ∧ M ),
meaning that we do not count the spikes that may occur after t but before s. This can be generalized to more general point processes. The main simulation loop is presented in Algorithm 1

Algorithm 1 Classical point process simulation algorithm

⊲ With [t 0 , t 1 ] the interval of simulation 1: Initialize the family of points P = ∅ ⊲ Each point is a time T with a mark, j T , which is the neuron on which T appears 2: Initialize t ← t 0 3: repeat 4:

for each neuron i ∈ I do

5:

Draw independently an exponential variable E i with parameter 1

6:

Apply the inverse transformation, that is, find T i such that

T i t φ abs i (s, t)ds = E i . 7:
end for 8:

Compute the time T of the next spike of the system after t, and the neuron where the spike occurs by T ← min i∈I T i , with j T ← arg min i∈I T i 9:

if T ≤ t 1 then 10:

append T with mark j T to P 11:

end if

12: t ← T 13: until t > t 1
Note that the quantity φ abs i (s, t) can be also seen as the hazard rate of the next potential point T

(1) i after t. It is a discrete event approach with the state corresponding to the function φ abs i (., t). Ogata [START_REF] Ogata | On Lewis Simulation Method for Point Processes[END_REF], inspired by Lewis' algorithm [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes[END_REF], added a thinning (also called rejection) step on top of this procedure because the integral

T (1) i t φ abs i (s, t)
ds can be very difficult to compute. To do so (and simplifying a bit), assume that φ i (t|F int t-) is upper bounded by a deterministic constant M . This means that the point process has always less points than a homogeneous Poisson process with intensity M . Therefore Steps 5-6 of Algorithm 1 can be replaced by the generation of an exponential of parameter M , E ′ i and deciding whether we accept or reject the point with probability φ abs i (t + E ′ i , t)/M . There are a lot of variants of this procedure: Ogata's original one uses actually the fact that the minimum of exponential variables, is still an exponential variable. Therefore one can propose a next point for the whole system, then accept it for the whole system and then decide on which neuron of the network the event is actually appearing. More details on the multivariate Ogata's algorithm can be found in [START_REF] Mascart | Centralized and distributed simulations of point processes using local independence graphs: A computational complexity analysis[END_REF].

As we see here, Ogata's algorithm is very general but clearly needs to simulate the whole system to simulate only one neuron. Moreover starting at time t 0 , it does not go backward and therefore cannot simulate a Hawkes process in stationary regime. There has been specific algorithms based on clusters representation that aim at perfectly simulate particular univariate Hawkes processes [START_REF] Møller | Perfect simulation of Hawkes processes[END_REF]. The algorithm that we propose here, will also overcome this flaw.

Kalikow decomposition

Kalikow decomposition relies on the concept of neighborhood, denoted by v, which are picked at random and which gives the portion of time and neuron subsets that we need to look at, to move forward. Typically, for a positive constant A, such a v can be:

• {(i, [-A, 0))}, meaning we are interested only by the spikes of neuron i in the window [-A, 0);

• {(i, [-2A, 0)), (j, [-2A, -A))}, that is, we need the spikes of neuron i in the window [-2A, 0)
and the spikes of neuron j in the window [-2A, -A); • the emptyset ∅, meaning that we do not need to look at anything to pursue.

We need to also define l(v) the total time length of the neighborhood v whatever the neuron is. For instance, in the first case, we find l(v) = A, in the second l(v) = 3A and in the third l(v) = 0.

We are only interested by stationary processes, for which the conditional intensity, φ i (t | F int t -), only depends on the intrinsic distance between the previous points and the time t and not on the precise value of t per se. In this sense the rule to compute the intensity may be only defined at time 0 and then shifted by t to have the conditional intensity at time t. In the same way, the timeline of a neighborhood v is defined as a subset of R * -so that information contained in the neighborhood is included in F int 0-, and v can be shifted (meaning its timeline is shifted) at position t if need be. We assume that I the set of neurons is countable and that we have a countable set of possibilities for the neighborhoods V.

Then, we say that the process admits a Kalikow decomposition with bound M and neighborhood family V, if for any neuron i ∈ I, for all v ∈ V there exists a non negative M -bounded quantity φ v i , which is F int 0-measurable and whose value only depends on the points appearing in the neighborhood v, and a probability density function λ i (.) such that

φ i (0 | F int 0 -) = λ i (∅)φ ∅ i + v∈V,v =∅ λ i (v) × φ v i (3.1) with λ i (∅) + v∈V,v =∅ λ i (v) = 1.
Note that because of the stationarity assumptions, the rule to compute the φ v i 's can be shifted at time t, which leads to a predictable function that we call φ vt i (t) which only depends on what is inside v t , which is the neighborhood v shifted by t. Note also that φ ∅ i , because it depends on what happens in an empty neighborhood, is a pure constant.

The interpretation of (3.1) is tricky and is not as straightforward as in the discrete case (see [START_REF] Ost | Sparse space-time models: concentration inequalities and Lasso[END_REF]). The best way to understand it is to give the theoretical algorithm for simulating the next event on neuron i after time t (cf. Algorithm 2).

Algorithm 2 Kalikow theoretical simulation algorithm

⊲ With [t 0 , t 1 ] the interval of simulation for neuron i 1: Initialize the family of points P = ∅ ⊲ NB: since we are only interested by points on neuron i, j T = i is a useless mark here. Draw an exponential variable E with parameter M , and compute T = t + E.

5:

Pick a random neighborhood according to the distribution λ i (.) given in the Kalikow decomposition and shift the neighborhood at time T : this is V T .

6:

Draw X T a Bernoulli variable with parameter φ

V T i (T ) M 7:
if X T = 1 and T ≤ t 1 then 8:

append T to P 9:

end if

10: t ← T 11: until t > t 1
This Algorithm is close to Algorithm 1 but adds a neighborhood choice (Step 5) with a thinning step (Steps 6-9).

In Appendix A, we prove that this algorithm indeed provides a point process with an intensity given by (3.1) shifted at time t.

The previous algorithm cannot be put into practice because the computation of φ VT i depends on the points in V T , that are not known at this stage. That is why the efficient algorithm that we propose in the next section goes backward in time before moving forward.

Backward Forward algorithm

Let us now describe the complete Backward Forward algorithm (cf. Algorithm 3). Note that to do so, the set of points P is not reduced, as in the two previous algorithms, to the set of points that we want to simulate but this contains all the points that need to be simulated to perform the task.

Algorithm 3 Backward Forward Algorithm

⊲ With [t 0 , t 1 ] the interval of simulation for neuron i ∈ I 1: Initialize the family V of non empty neighborhoods with {(i, [t 0 , t 1 ])} 2: Initialize the family of points P = ∅ ⊲ Each point is a time T with 3 marks: j T is the neuron on which T appears, V T for the choice of neighborhood, X T for the thinning step (accepted/rejected) 3: Draw E an exponential variable with parameter M 4: Schedule Tnext = t 0 + E 5: while Tnext < t 1 do

6:

Append to P, the point Tnext, with 3 marks: j Tnext = i, V Tnext = n.a. and X Tnext = n.a. (n.a. stand for not assigned yet)

7:

while There are points T in P with V T = n.a. do

8:

for each point T in P with V T = n.a. do

9:

Update V T by drawing V T according to λ j T shifted at time T . 10:

if V T = ∅ then 11:
Find the portion of time/neurons in V T which does not intersect the existing non empty neighborhoods in V 12:

Simulate on it a Poisson process with rate M 13:

Append the simulated points, T ′ , if any, to P with their neuron j T ′ and with V T ′ = X T ′ =n.a Sort the T 's in P with X T =n.a. in increasing order

19:

for each of them starting with the most backward do Tnext ← Tnext + E ′ 24: end while 25: The desired points are the points in P with marks j T = i, X T = 1 and that appear in [t 0 , t 1 ]

⊲ It is possible that the algorithm generated points before t 0 and they have to be removed Init.

Init. Marks

Backward

Forward

At the difference with Algorithm 2, the main point is that in the backward part we pick at random all the points that may influence the thinning step. The fact that this loop ends comes from the following Proposition.

Proposition 1. If sup i∈I v∈V λ i (v)l(v)M < 1. (4.1)
then the backward part of Algorithm 3 ends almost surely in finite time.

The proof is postponed to Appendix B. It is based on branching process arguments. Basically if in Steps 8-16, we produce in average less than one point, either because we picked the empty set in V T or because the simulation of the Poisson process ended up with a small amount of points, eventually none, then the loop ends almost surely because there is an extinction of the corresponding branching process.

In the backward part, one of the most delicate part consists in being sure that we add new points only if we have not visited this portion of time/neurons before (see ). If we do not make this verification, we may not have the same past depending on the neuron we are looking at and the procedure would not simulate the process we want.

In the forward part, because the backward algorithm stopped just before, we are first sure to have assess all V T 's. Since φ Vt j (t) is F int t-measurable, for all t, φ VT jT (T ) only depends on the points in P with mark X T = 1 inside V T . The problem in Algorithm 2, phrased differently, is that we do not know the marks X T of the previous points when we have to compute φ VT jT (T ). But in the forward part of Algorithm 3, we are sure that the most backward point for which the thinning (X T =n.a.) has not taken place, satisfies

• either V T = ∅ • or V T = ∅ but
either there are no simulated points in the corresponding V T or the points there come from previous rounds of the loop (Step 5). Therefore their marks X T have been assigned.

Therefore, with the Backward Forward algorithm, and at the difference to Algorithm 2, we take the points in an order for which we are sure that we know the previous needed marks. Figure 1 describes an example to go step by step through Algorithm 3. The backward steps determine all the points that may influence the acceptation/rejection of point T next . Notice that whereas usual activity tracking algorithms for point processes [START_REF] Mascart | Centralized and distributed simulations of point processes using local independence graphs: A computational complexity analysis[END_REF] automatically detect the active children (influencees), activity tracking in the backward steps detect the parents (influencers). The forward steps finally select the points.

Illustration

To illustrate in practice the algorithm, we have simulated a Hawkes process as given in (2.1). Indeed such a process has a Kalikow decomposition (3.1) with bound M and neighborhood family V constituted of the v's of the form v = {(j, [-A, 0))} for some neuron j in I. To do that, we need the following choices:

λ i (∅) = 1 - j∈I w ij and φ ∅ i = ν i λ i (∅)
and for v of the form v = {(j, [-A, 0))} for some neuron j in I,

λ i (v) = w ij and φ v i = Nb j [-A,0) ∧ M.
We have taken I = Z 2 and the w ij proportional to a discretized centred symmetric bivariate Gaussian distribution of standard deviation σ. More precisely, once λ i (∅) = λ ∅ fixed, picking according to λ i consists in

• choosing whether V is empty or not with probability λ ∅ • if V = ∅, choosing V = {(j, [-A, 0))} with j -i = round(W ) and W obeys a bivariate N (0, σ 2 ).
In all the sequel, the simulation is made for neuron i = (0, 0) with t 0 = 0, t 1 = 100 (see Algorithm 3). The parameters M, λ ∅ and σ vary. The parameters ν i = ν and A are fixed accordingly by

ν = 0.9M λ ∅ and A = 0.9M -1 (1 -λ ∅ ) -1 ,
to ensure that φ ∅ i < M and (4.1), which amounts here to (1λ ∅ )AM < 1. On Figure 2(a), with M = 2, σ = 1 and λ ∅ small, we see the overall spread of the algorithm around the neuron to simulate (here (0, 0)). Because we chose a Gaussian variable with small variance for the λ i 's, the spread is weak and the neurons very close to the neuron to simulate are requested a lot of time at Steps 9-11 of the algorithm. This is also where the algorithm spent the biggest amount of time to simulate Poisson processes. Note also that roughly to simulate 80 points, we need to simulate 10 times more points globally in the infinite network. Remark also on Figure 2(b), the avalanche phenomenon, typical of Hawkes processes: for instance the small cluster of black points on neuron with label 0 (i.e. (0,0)) around time 22, is likely to be due to an excitation coming for the spikes generated (and accepted) on neuron labeled 8 and self excitation. The beauty of the algorithm is that we do not need to have the whole time line of neuron 8 to trigger neuron 0, but only the small blue pieces: we just request them at random, depending on the Kalikow decomposition.

On Figure 3, we can first observe that when the parameter which governs the range of λ i 's increase, the global spread of the algorithm increase. In particular, comparing the top left of Figure 3 to Figure 2 where the only parameter that changes is σ, we see that the algorithm is going much further away and simulates much more points for a sensible equal number of points to generate (and accept) on neuron (0,0). Moreover we can observe that • From left to right, by increasing λ ∅ , it is more likely to pick an empty neighborhood and as a consequence, the spread of the algorithm is smaller. By increasing ν = 0.9M λ ∅ , this also increases the total number of points produced on neuron (0;0). • From top to bottom, by increasing M , there are more points which are simulated in the Poisson processes (Step 12 of Algorithm 3) and there is also a stronger interaction (we do not truncate that much the number of points in φ v ). Therefore, the spread becomes larger and more uniform too, because there are globally more points that are making requests. Moreover, by having a basic rate M which is 10 times bigger, we have to simulate roughly 10 times more points. 8 ). Note that since the simulation is on [0, 100] the time spent at position (0, 0) is at least 100. On (a), the summary for one simulation with below the plot, number of points accepted at neuron (0, 0) and total number of points that have been simulated. Also annotated on (a), with labels between 0 and 8, the 9 neurons for which the same simulation in [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF]40] is represented in more details on (b). More precisely on (b), in abscissa is time and the neuron labels in ordinate. A plain dot represents an accepted point, by the thinning step (Step 20 of Algorithm 3), and an empty round, a rejected point. The blue pieces of line represent the non empty neighborhoods that are in V.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 1 2 3 4 5 6
• • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • (
• • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -20 0 20 -20 0 20
Accepted at (0,0): 86 Produced : 2001

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -30

Conclusion

We derived a new algorithm for simulating the behavior of one neuron embedded in an infinite network. This is possible thanks to the Kalikow decomposition which allows picking at random the influencing neurons. As seen in the last section, it is computationnally tractable in practice to simulate open systems in the physical sense. A question that remains open for future work is whether we can prove that such a decomposition exists for a wide variety of processes, as it has been shown in discrete time (see [START_REF] Galves | Infinite Systems of Interacting Chains with Memory of Variable Length A Stochastic Model for Biological Neural Nets[END_REF][START_REF] Galves | Modeling networks of spiking neurons as interacting processes with memory of variable length[END_REF][START_REF] Ost | Sparse space-time models: concentration inequalities and Lasso[END_REF]).

thatis F N (1) t-measurable and therefore F N t-measurable . We have, for any i ∈ I,

E   ∞ 0 C t dN i t (1)   = ∞ n=1 E C T i n 1 X i n =1
Note that by Theorem T35 at Appendix A1 of [START_REF] Brémaud | Point Processes and Queues: Martingale Dynamics[END_REF], any point T should be understood as a stopping time, and that by Theorem T30 at Appendix A2 of [START_REF] Brémaud | Point Processes and Queues: Martingale Dynamics[END_REF],

F N T -= j σ{T j m , X j m such that T j m < T } So E   ∞ 0 C t dN i t (1)   = ∞ n=1 E C T i n E(1 X i n =1 |F N T i n -, V i n ) = ∞ n=1 E C T i n ϕ V i n i (T i n ) M .
Let us now integrate with respect to the choice V i n , which is independent of anything else.

E   ∞ 0 C t dN i t (1)   = ∞ n=1 E   CT i n λ i (∅)ϕ ∅ i + v∈V,v =∅ λ i (v) × ϕ v i (T i n ) M    = E   ∞ 0 C t ϕ i (
  = E   ∞ 0 C t ϕ i (t)dt   ,
which ends the proof.

Therefore, if we denote the total number of sites in Ck by Z (k) , we have

E(Z (k+1) | Ck ) ≤ Z (k) sup i∈I ζ i .
One can then conclude by recursion that,

E(Z (k) ) ≤ (sup i∈I ζ i ) k < 1.
The last inequality use the sparsity neighborhood assumption. Then we deduce that, the mean number of children in each generation goes to 0 as k tends to infinity. So by using classical branching techniques in [START_REF] Méléard | Aléatoire: introduction à la théorie et au calcul des probabilités[END_REF], we conclude that the tree C will go extinct almost surely. This also implies that, the backward steps end a.s.

2 :

 2 Initialize t ← t 0 3: repeat 4:

20 :

 20 Draw X T as a Bernoulli variable with parameter φ Draw E ′ another exponential variable with parameter M 23:

Fig 1 : 3 ,

 13 Fig 1: Main flow example for Algorithm 3, with backward steps in green (cf. Algorithm 3, Steps 7-17) and forward steps in purple (cf. Algorithm 3, Steps 18-25). Following arrow numbers: (1) The next point T next = t + E (cf. Algorithm 3, Step 4) is scheduled, (2) The neighborhood V Tnext is selected in the first backward step, a first generation of three points (a, b on neuron k and c on neuron ℓ) is drawn (cf. Algorithm 3, Step 9), thanks to a Poisson process, (cf. Algorithm 3, Steps 11-12) and appended to P (cf. Algorithm 3, Step 13), (3) at the second generation, a non empty neighborhood is found, i.e. V b = ∅ (cf. Algorithm 3, Steps 9-1), but the Poisson process simulation does not give any point in it (cf. Algorithm 3, Step 12), (4) at the second generation, the neighborhood V a is picked, it is not empty and overlap the neighborhood of the first generation (cf. Algorithm 3, Steps 9-11): therefore there is no new simulation in the overlap (c is kept and belongs to V b as well as V a ) but there is a new simulation thanks to a Poisson process outside of the overlap leading to a new point d (cf. Algorithm 3, Step 12)(5) at the second generation, for point c, one pick the empty neighborhood, i.e. V c = ∅ (cf. Algorithm 3, Step 9) and therefore we do not simulate any Poisson process, (6) at third generation, similarly no point and no interval are generated, i.e. V d = ∅ (cf. Algorithm 3, Step 9). This is the end of the backward steps and the beginning of the forward ones, (7) the point d is not selected, acceptation/selection taking place with probability φ ∅ ℓ M (cf. Algorithm 3, Step 20), (8) the point c is accepted, here again with probability φ ∅ ℓ M (cf. Algorithm 3, Step 20), (9) the point b is not selected, acceptation taking place, here, with probability φ V b k (b) M (cf. Algorithm 3, Step 20), (10) the point a is selected, acceptation taking place, here, with probability φ Va k (a) M (cf. Algorithm 3, Step 20), (11) The neighborhood of neuron i contains two points, one on neuron k and one on neuron ℓ and one selects T next with probability φ V T next i

to 9 Fig 2 :

 92 Fig 2:Simulation for M = 2, σ = 1, λ ∅ = 0.25. For each neuron in Z 2 , that have been requested in Steps 9:11, except the neuron of interest (0, 0), have been counted the total number of requests, that is the number of time a V T pointed towards this neuron (Steps 9 and 11) and the total time spent at this precise neuron simulating a homogeneous Poisson process (Step 12). Note that since the simulation is on [0, 100] the time spent at position (0, 0) is at least 100. On (a), the summary for one simulation with below the plot, number of points accepted at neuron (0, 0) and total number of points that have been simulated. Also annotated on (a), with labels between 0 and 8, the 9 neurons for which the same simulation in[START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF]40] is represented in more details on (b). More precisely on (b), in abscissa is time and the neuron labels in ordinate. A plain dot represents an accepted point, by the thinning step (Step 20 of Algorithm 3), and an empty round, a rejected point. The blue pieces of line represent the non empty neighborhoods that are in V.

Fig 3 :

 3 Fig 3: Simulation for 4 other sets of parameters, all of them with σ = 3. Summaries as explained in Figure 2. On top, M = 2; on bottom, M = 20. On the left part, λ ∅ = 0.25, on the right part, λ ∅ = 0.5.

  Since Π i is a Poisson process with respect to (F N t ) t with intensity M , and since C t

				
			t) M	dΠ i (t)  .
	measurable, we finally have that	ϕi(t) M	is F N t-
			∞
	E		C t dN i t (1)
			0

Acknowledgements

This work was supported by the French government, through the UCA Jedi Investissements d'Avenir managed by the National Research Agency (ANR-15-IDEX-01) and by the interdisciplinary Insti-tute for Modeling in Neuroscience and Cognition (NeuroMod) of the Université Côte d'Azur. The authors would like to thank Professor E.Löcherbach from Paris 1 for great discussions about Kalikow decomposition and Forward Backward Algorithm.

Accepted at (0,0): 154 Produced : 1204

Accepted at (0,0): 969 Produced : 15378 To prove that Algorithm 2 returns the desired processes, let us use some additional and more mathematical notation. Note that all the points simulated on neuron i before being accepted or not can be seen as coming from a common Poisson process of intensity M , denoted Π i . For any i ∈ I, we denote the arrival times of Π i , (T i n ) n∈Z , with T i 1 being the first positive time. As in Step 6 of Algorithm 2, we attach to each point of Π i a stochastic mark X given by,

Let us also define V i n the neighborhood choice of T i n picked at random and independently of anything else according to λ i and shifted at time T i n . In addition, for any i ∈ I, define N i = (T i n , X i n ) n∈Z an E-marked point process with E = {0; 1}. In particular, following the notation in Chapter VIII of [START_REF] Brémaud | Point Processes and Queues: Martingale Dynamics[END_REF], for any i ∈ I, let

Moreover note that (N i t (1)) t∈R is the counting process associated to the point processs P simulated by Algorithm 2. Let us denote by ϕ i (t), the formula given by (3.1) and shifted at time t. Note that since the φ v i 's are Proof. We do the proof for the backward part, starting with T = T next as the next point after t 0 (Step 4 of Algorithm 3), the proof being similar for the other T next generated at Step 23. We construct a tree with root (i, T ). For each point (j T ′ , T ′ ) in the tree, the points which are simulated in V T ′ (Step 12 of Algorithm 3) define the children of (j T ′ , T ′ ) in the tree. This forms the tree T .

Let us now build a tree C with root (i, T ) (that includes the previous tree) by mimicking the previous procedure in the backward part, except that we simulate on the whole neighborhood even if it has a part that intersects with previous neighborhoods (if they exist) (Step 11-12 of Algorithm 3). By doing so, we make the number of children at each node independent of anything else.

If the tree C goes extinct then so does the tree T and the backward part of the algorithm terminates.

But if one only counts the number of children in the tree C, we have a marked branching process whose reproduction distribution for the mark i is given by • no children with probability λ i (∅)

• Poissonian number of children with parameter l(v)M if v is the chosen neighborhood with probability λ i (v)

This gives that the average number of children issued from a node with the mark i is

If we denote Ck as the collection of points in the tree C at generation k, and by K T ′ the set of points generated independently as a Poisson process of rate M inside V T ′ , we see recursively that Ck+1 =