Emilie Coupechoux
email: emilie.coupechoux@ens.fr

Marc Lelarge
email: marc.lelarge@ens.fr

Impact of Clustering on Diffusions and Contagions in Random Networks

Keywords: Contagion threshold, diffusion, Random graphs, clustering

Motivated by the analysis of social networks, we study a model of network that has both a tunable degree distribution and a tunable clustering coefficient. We compute the asymptotic (as the size of the population tends to infinity) for the number of acquaintances and the clustering for this model. We analyze a contagion model with threshold effects and obtain conditions for the existence of a large cascade. We also analyze a diffusion process with a given probability of contagion. In both cases, we characterize conditions under which a global cascade is possible.

I. INTRODUCTION

Most of the epidemic models [START_REF] Newman | Spread of epidemic disease on networks[END_REF], [START_REF] Pastor-Satorras | Evolution and Structure of the Internet: A Statistical Physics Approach[END_REF] consider a transmission mechanism which is independent of the local condition faced by the agents concerned. There is now a vast literature on epidemics in complex networks (see [START_REF] Newman | Properties of highly clustered networks[END_REF] for a review) and there is now a good understanding of the impact of the topology on the spread of an epidemic. But if there is a factor of persuasion or coordination involved, relative considerations tend to be important in understanding whether some new behavior or belief is adopted [START_REF] Vega-Redondo | Complex social networks[END_REF].

In social contexts, the diffusion of information and behavior often exhibits features that do not match well those of the SIR or SIS model [START_REF] Vega-Redondo | Complex social networks[END_REF]. In the classical (SI) diffusion model, an individual is influenced by each of her neighbors independently. For the spread of a new technology in the network, a rather appropriate model is the case where each individual adopts the technology as soon as enough of her neighbors have already adopted it: this corresponds to the basic game-theoretic contagion model proposed by Morris [10]. Consider a graph G in which the nodes are the individuals in the population and there is an edge (i, j) if i and j can interact with each other. Each node has a choice between two possible behaviors labeled A and B. On each edge (i, j), there is an incentive for i and j to have their behaviors match, which is modeled as the following coordination game parametrized by a real number q ∈ (0, 1): if i and j choose A (resp. B), they each receive a payoff of q (resp. (1q)); if they choose opposite strategies, then they receive a payoff of 0. Then the total payoff of a player is the sum of the payoffs with each of her neighbors. If the degree of node i is d i and S B i is the number of its neighbors playing B, then the payoff to i from choosing A is q(d i -S B i ) while the payoff from choosing B is (1q)S B i . Hence, in best response update, i should adopt B if S B i > qd i and A if S B i ≤ qd i . A number of qualitative insights can be derived from a diffusion model even at this level of simplicity. Specifically, consider a network where all nodes initially play A. If a small number of nodes are forced to adopt strategy B (the seed) and we apply best-response updates to other nodes in the network, then these nodes will be repeatedly applying the following rule: switch to B if enough of your neighbors have already adopted B. There can be a cascading sequence of nodes switching to B such that a network-wide equilibrium is reached in the limit.

Large complex networks such as social contact structures, the internet and various types of collaboration networks have received a lot of attention during the last few years; [START_REF] Newman | The Structure and Dynamics of Networks[END_REF] and the references therein. As for social networks, one of their most striking features is that they are highly clustered, meaning that there is a large number of triangles and other short cycles [START_REF] Newman | Properties of highly clustered networks[END_REF]. This is a consequence of the fact that friendship circles are typically strongly overlapping so that many of our friends are also friends of each other. A model (inspired from [START_REF] Trapman | On analytical approaches to epidemics on networks[END_REF]) that captures this in a natural way will be described in Section II. Roughly, the idea of the model is to 'add' clustering to a standard configuration model by replacing some vertices by cliques. By choosing the fraction of vertices replaced, this leads to a graph where the amount of clustering can be tuned by adjusting the parameters of the model. As we will show this model generalizes the standard configuration model to incorporate clustering and it is still possible to derive rigorously exact formulas for the analysis of contagions and diffusions on these networks. The model has the advantage to allow any arbitrary degree distribution: in particular, it can be applied to scale-free networks that have a power law degree distribution.

An important goal of network modeling is to investigate how the structure of the network affects the behavior of various types of dynamic processes on the network. The aim of this paper is to give a rigorous analysis of how clustering in a network affects the spread of an epidemic (we study both gametheoretical contagion and classical diffusion models). For the Reed-Frost epidemic, [START_REF] Newman | Properties of highly clustered networks[END_REF] studies the effect of clustering for a different model than ours and by heuristic means. Calculations indicate that the epidemic threshold should decrease as the clustering increases. This result has been recently rigorously proven in [START_REF] Britton | Epidemics on random graphs with tunable clustering[END_REF]. The analysis of a contagion process with threshold effect has not been done for their model. Another model of random graphs with positive clustering (different than the one we consider here) is introduced in [START_REF] Newman | Random graphs with clustering[END_REF], and heuristic results on both diffusion and contagion models are present in [START_REF] Hackett | Cascades on a class of clustered random networks[END_REF], for this model. Up to our knowledge, there is no rigorous analysis for the contagion model on a random graph with clustering.

Our main contributions are

• the study of a tractable random graph model with tunable clustering (inspired from [START_REF] Trapman | On analytical approaches to epidemics on networks[END_REF]). In Section II, we introduce it and compute the asymptotic degree distribution of the graph and the clustering coefficient, • the analysis of a contagion process with threshold effect.

In Section III, we derive the contagion threshold for our random graph model with clustering extending recent results of [START_REF] Lelarge | Diffusion and cascading behavior in random networks[END_REF], and • the analysis of a diffusion process with given probability of infection. In Section IV, we derive the minimal value for the probability of contagion in our random graph model with clustering such that a global diffusion is possible. This proves a heuristic result of [START_REF] Gleeson | Analytical results for bond percolation and k-core sizes on clustered networks[END_REF]. In the following, we consider asymptotics as n → ∞, and we denote by → p the convergence in probability as n → ∞. The abbreviation 'whp' ("with high probability") means with probability tending to 1 as n → ∞, and we use the notation o p (n) in a standard way: X = o p (n) means that, for every ε > 0, P(X > εn) → 0 as n → ∞.

II. RANDOM GRAPH MODEL AND ITS BASIC PROPERTIES

We first present the model for the social graph, then its asymptotic degree distribution, and finally its clustering coefficient.

A. Model

Let n ∈ N and let d = (d

(n) i ) n i=1 = (d i ) n
1 be a sequence of non-negative integers such that i d i is even. Let G (n, d) be a graph chosen uniformly at random among all graphs with n vertices and degree sequence d (assuming there exists such a graph) [START_REF] Bollobás | Random graphs[END_REF].

We will let n → ∞ and assume that we are given d satisfying the following regularity conditions, see [START_REF] Molloy | A critical point for random graphs with a given degree sequence[END_REF]:

Condition 1. For each n, d = (d i ) n
1 is a sequence of nonnegative integers such that i d i is even. We assume that there exists a probability distribution p = (p r ) ∞ r=0 (independent of n) such that: The model of random graphs G (n, d) has the advantage to allow to handle arbitrary degree distributions. However, these graphs are 'locally tree-like', i.e. they contain no short loops in their structure. We now show that it is possible to generalize random graphs to incorporate clustering in simple fashion and still to derive rigorously exact formulas for diffusions and contagions. The model of random graphs is based on the model G (n, d), but we 'add' clustering. The idea is to replace some vertices by a clique of size the degree in the original graph, i.e. a vertex of degree r in the original graph G (n, d) is replaced by r vertices with all the r(r -1)/2 edges between them and each of them is connected to exactly one of the neighbors of the vertex in the original graph G (n, d) as illustrated on the figure:

(i) n r /n = |{i : d i = r}| /n → p r as n → ∞, for all r ≥ 0; (ii) λ := r rp r ∈ (0, ∞); (iii) i d 3 i = O(n). If D n is
In order to be able to tune the clustering coefficient in the graph, we will not replace each vertex by a clique but rather do a probabilistic choice whether to replace a vertex or not. The resulting random graph will be denoted by G (n, d, γ), where γ = (γ r ) ∞ r=0 is a sequence such that: for all r ≥ 0, γ r ∈ [0, 1] represents the probability that a vertex of degree r in G (n, d) is replaced by a clique of size r in the new model G (n, d, γ). More precisely, for each vertex i ∈ {1, ..., n}, let X(i) be a Bernoulli variable with parameter γ di (all Bernoulli variables being independent). We construct the random graph G (n, d, γ) the following way: start from G (n, d) and, for each i ∈ {1, ..., n}, if X(i) = 1, replace i by a clique of size d i where each vertex of the clique has exactly one neighbor outside the clique being a neighbor of i in the original graph. Note in particular that if we choose γ r = 0 for all r ≥ 0, then G (n, d, γ) = G (n, d), whereas for γ r = 1 for all r ≥ 0, all vertices in G (n, d) have been replaced by cliques.

B. Degree distribution in G (n, d, γ)

As we will see in the next subsection, the procedure described above introduces clustering at soon as γ r > 0 for some r. It also modifies the degree distribution in the graph and we derive the new degree distribution here. Let

i ∈ {1, ..., n}. If X(i) = 1, vertex i in G (n, d) is replaced in G (n, d, γ
) by a clique of d i vertices, all having degree d i (indeed, each vertex of the clique has d i -1 edges linked with the other vertices of the clique, and one 'external' edge). So each vertex of degree r can either be replaced (with probability γ r ) by r vertices of degree r, or stays as a single vertex of degree r (with probability 1γ r ). The following proposition gives the asymptotic degree distribution in G (n, d, γ). Proposition 2. We consider the model G (n, d, γ) for a sequence d satisfying Condition 1 with probability distribution p = (p r ) ∞ r=0 , and clustering parameter γ = (γ r ) ∞ r=0 . For all r ≥ 0, let ñr be the number of vertices with degree r in G (n, d, γ), and let ñ = r ñr be the total number of vertices in G (n, d, γ). Then we have, as n → ∞:

ñ n p -→ γ := d≥0 [dγ d + (1 -γ d )] p d
and, for all r ≥ 0, the proportion of vertices with degree r in G (n, d, γ) has the following limit, as n → ∞: Note that the number of vertices with degree d that are not replaced by a clique is n d -B d , so we can express the number ñ of vertices in G (n, d, γ) the following way:

ñr ñ p -→ pr := [rγ r + (1 -γ r )] p r γ . Proof: Let d ≥ 0,
ñ n = 1 n d dB d + (n d -B d ) p -→ d [dγ d + (1 -γ d )] p d = γ
which follows from the previous limits, and the uniform integrability of the random variables

D n (see Condition 1- (iii)).
In particular, this shows that the total number of vertices with degree r (r ≥ 0) in G (n, d, γ) is rB r + (n r -B r ), so the proportion of vertices with degree r in G (n, d, γ) is:

rB r + (n r -B r ) ñ p -→ [rγ r + (1 -γ r )] p r γ = pr ,
which concludes the proof.

In other words, if Dn is the degree of a vertex chosen uniformly at random in G (n, d, γ), then Proposition 2 implies that Dn d -→ D, where D is a random variable with distribution (p r ) r≥0 .

In our definition of G (n, d, γ), each vertex of degree 0 in G (n, d) is removed from the graph with probability γ 0 , and kept with probability 1γ 0 . We could have considered the case where each vertex of degree 0 is kept with probability 1.

In that case, we have that p0 = p 0 and the following value γ: γ = d≥1 [dγ d + (1γ d )] p d + p 0 . To simplify notations, we keep the definition of the Proposition, but arguments are the same.

In the particular case where γ r = γ for all r, we have γ = γλ + 1γ and the mean degree of D is then

λ = E[ D] = γE[D 2 ] + (1 -γ)λ γλ + 1 -γ ,
which is a non-decreasing function of γ.

C. Clustering coefficient

The local clustering coefficient

C (n) v
of a vertex v in a graph quantifies how close the vertex and its neighbors are to being a clique (complete graph) [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF]. C (n) v is defined to be the fraction of pairs of neighbors of v that are neighbors also of each other. More formally, let N v be the set of neighbors of v (its cardinality |N v | = d v is the degree of v), and let P v be the number of pairs {w, w ′ } ⊂ N v , w = w ′ , such that w and w ′ share an edge together. The total number of possible pairs is d v (d v -1)/2, so we define the local clustering coefficient of v as

C (n) v = P v • 2/[d v (d v -1)
]. The biased clustering coefficient for the whole network C (n) is defined as the average of the local clustering coefficient for each vertex:

C (n) = v C (n) v /ñ,
C (n) p -→ C := r≥3 p r γ r γ (r -2) ,
where γ is defined in Proposition 2.

The proof is given at the end of the subsection (together with the proof of Proposition 4).

In the particular case where γ r = γ ≥ 0 for all r ≥ 0, we get for the asymptotic biased clustering coefficient: of the graph by:

C = λ -2 + p 1 + 2p 0 λ -1 + 1 γ . If γ = 0,
C (n) 2 = 2 • v P v v d v (d v -1)
.

It is the mean probability that three given vertices constitute a triangle conditional on that two of the three possible edges between them exist, and it corresponds to the notion of clustering studied in [START_REF] Britton | Epidemics on random graphs with tunable clustering[END_REF].

C (n) 2 p -→ C 2 := r≥2 r(r -1)(r -2)γ r p r r≥2 ((r -1)γ r + 1)r(r -1)p r .

Proofs of Propositions 3 and 4:

We have the following result, that follows from the fact that G (n, d) converges locally to a tree, when n → ∞:

Lemma 5. Let C (n) be the biased clustering coefficient in G (n, d).
Then we have:

C (n) p -→ 0.
The same result holds for the clustering coefficient

C (n) 2 of the graph G (n, d).
We say that a vertex in G (n, d, γ) has parent i ∈ {1, ..., n} it belongs to a clique that replaces the vertex i of G (n, d) (when

X(i) = 1) or if it is i (when X(i) = 0).
We first consider a vertex v in G (n, d, γ) whose parent i is such that

X(i) = 1 (vertex i of G (n, d) is replaced by a clique K in G (n, d, γ)).
In this case we can directly compute the local clustering coefficient C (n) v . Indeed, vertex v has d i -1 neighbors inside K, that are all linked together (which gives

(di-1)(di-2) 2
edges in total), and one neighbor v ′ outside K, which is not linked to the other neighbors of v (if it were the case, there would be several edges between i and the parent j of v ′ , which is not the case in the simple graph G (n, d)). Hence

P v = (d i -1)(d i -2) 2 and C (n) v = 2P v d i (d i -1) = d i -2 d i , provided that d i ≥ 2. If d i ∈ {0, 1}, then C (n) v
= 0. We first prove Proposition 3. Since there are d i such vertices inside a clique, the contribution of clique K in the total clustering

C (n) = v C (n) v /ñ is equal to d i C (n) v /ñ = (d i -2)/ñ.
This leads to the following:

ñ n C (n) = 1 n d≥2 (d -2)B d + 1 n i:X(i)=0 C (n) i
where B d is the number of vertices with degree d that are replaced by a clique, as in the proof of Proposition 2. Using that B d /n → p γ d p d , and that i:X(i)=0 C

(n) i /n → p 0 (as a consequence of Lemma 5), we obtain:

ñ n C (n) p -→ d≥3 (d -2) γ d p d . Proposition 3 follows, applying Proposi- tion 2.
The end of the proof for Proposition 4 is similar, and follows from the fact that:

2 v P v /n → p d d(d -1)(d -2)γ d p d , v d v (d v -1)/n → p d d(d -1)[dγ d + (1 -γ d )]p d .
In the next paragraph, we use the second definition of clustering coefficient, but results are true for both.

D. Link between clustering coefficient and mean degree in the graph G (n, d, γ)

Both asymptotic clustering coefficient C 2 (Proposition 4) and asymptotic degree distribution (p r ) r≥0 (Proposition 2) depend on the following parameters:

• the asymptotic degree distribution p of the original graph G (n, d), • the clustering parameter γ. We will see how the clustering coefficient C 2 and the mean degree λ = r r pr vary when p and γ vary. We focus on the case when γ r = γ for all r ≥ 0, i.e. each vertex in G (n, d) is replaced by a clique with probability γ. In addition, we impose some conditions on the probability distribution p.

Indeed, in figures 1 and 2, we assume that the initial graph G (n, d) has power law degree distribution with exponential cutoff: there exists a power τ > 0 and a cutoff κ > 0 such that, for all r ≥ 1, p r = c(τ, κ) • r -τ e -r/κ , where c(τ, κ) = 1/( s s -τ e -s/κ ) is a normalizing constant. This cutoff κ allows Condition 1 to be satisfied for any power τ > 0: in figures, we take κ = 50. If we compute (using Proposition 2) the asymptotic degree distribution in G (n, d, γ) when each vertex is replaced by a clique (γ = 1), we obtain a power law distribution with parameter τ -1: pr ∝ r -(τ -1) e -r/κ . More generally, for any given value of γ, the degree distribution in G (n, d, γ) is a linear combination between a power law of parameter τ , and a power law of parameter τ -1.

Once we impose this form for the probability distribution p, we are left with two degrees of freedom: the power τ and the probability γ. In figure 1, we make these two parameters vary, and plot the correlation between the clustering coefficient C 2 , the mean degree λ and the probability γ.

Figure 2 represents several slices of figure 1, for different values of the clique probability γ. Increasing the mean degree in the graph G (n, d, γ) also increases the clustering. With our graph model, we are not able to reach a clustering coefficient of 1, especially for low values of the mean degree λ. Yet the maximal clustering coefficient we can obtain (for γ = 1) is greater than 0.9 as soon as the mean degree λ is greater than 2.5. Now we are interested in the evolution of the clustering coefficient C 2 with respect to γ, when the mean degree λ is fixed (which corresponds to other slices of figure 1). In order to keep λ fixed when γ varies, we need to adjust the value of the power τ : this changes both the probability distribution p and the asymptotic degree distribution p (whose mean is kept fixed).

For technical reasons, it will be more convenient to work with a Poisson distribution with mean λ > 0 for the probability distribution p: for all r ≥ 0, p r = e -λ λ r /r!, that is to say G (n, d) is distributed as an Erdős-Renyi graph. If γ = 1, the asymptotic degree distribution in G (n, d, γ) is a Poisson variable (with parameter λ) shifted by 1, so the mean degree is λ = λ + 1. In that case, the degrees of freedom are parameters λ and γ. In order to keep λ fixed when γ varies, we have to adjust the value of λ. If the mean degree λ in G (n, d, γ) is fixed and high enough, the clustering coefficient C 2 increases with γ (figure 3), but, for low values of λ, the clustering coefficient C 2 is not a non-decreasing function of γ (distributions p and p vary).

The advantage of this model is that it allows to consider any degree distribution, contrary to random intersection graphs [START_REF] Britton | Epidemics on random graphs with tunable clustering[END_REF] that are restricted to the Poisson distribution. Yet we are limited by the correlation between the clustering coefficient and the mean degree in the graph G (n, d, γ): in [START_REF] Britton | Epidemics on random graphs with tunable clustering[END_REF], the clustering coefficient can vary between 0 and 1, even for low values of the mean degree in the graph.

III. CONTAGION THRESHOLD FOR RANDOM GRAPHS WITH CLUSTERING

A. Contagion model

Motivated by the game-theoretic contagion model proposed by Morris [10] and described in the introduction, we now describe formally our model of contagion on any finite graph G. The progressive dynamic of the diffusion on the finite graph G operates as follows: some set of nodes S starts out being active; all other nodes are inactive. Time operates in discrete steps t = 1, 2, 3, . . . . At a given time t, any inactive node i becomes active if its number of active neighbors is at least ⌊qd i ⌋ + 1. This in turn may cause other nodes to become active. It is easy to see that the final set of active nodes (after n time steps if the network is of size n) only depends on the initial set S (and not on the order of the activations) and can be obtained as follows: set Y i = 1(i ∈ S) for all i. Then as long as there exists i such that j∼i Y j > qd i , set Y i = 1, where j ∼ i means that i and j share an edge in G. When this algorithm finishes, the final state of node i is represented by Y i : Y i = 1 if node i is active and Y i = 0 otherwise.

We see that the lower q is, the easier the diffusion spreads. In [10], the contagion threshold of a connected infinite network is defined as the maximum threshold q c at which a finite set of initial adopters can cause a complete cascade, i.e. the resulting cascade of adoptions of B eventually causes every node to switch from A to B. In this section, we restrict ourselves to the model where the initial adopters are forced to play B forever. In this case, the contagion is monotone and the number of nodes playing B is non-decreasing. We say that this case corresponds to the permanent adoption model: a player playing B will never play A again.

B. Phase transition for the contagion

We now compute the contagion threshold for a sequence of random networks. Since a random network is finite and not necessarily connected, we first need to adapt the definition of contagion threshold to our context as was done in [START_REF] Lelarge | Diffusion and cascading behavior in random networks[END_REF]. For a graph G = (V, E) and a parameter q, we consider the largest connected component of the induced subgraph in which we keep only vertices of degree strictly less than q -1 . We call the vertices in this component pivotal players: if only one pivotal player switches from A to B then the whole set of pivotal players will eventually switch to B in the permanent adoption model. For a player v ∈ V , we denote by C(v, q) the final number of players B in the permanent adoption model with parameter q, when the initial state consists of only v playing B, all other players playing A. Informally, we say that C(v, q) is the size of the cascade induced by player v. Theorem 6. Consider the random graph G (n, d, γ) for a sequence d satisfying Condition 1 with probability distribution p = (p r ) ∞ r=0 , and define

q c = sup    q : r<q -1 r(r -1)p r > r rp r    .
Let P(n) be the set of pivotal players of G (n, d, γ). (i) If q < q c , then there is a unique ξ ∈ (0, 1) such that

d<q -1 dp d (1 -ξ d-1 ) = λ(1 -ξ)
and we have:

| P(n) |/ñ p -→ d<q -1 [dγ d + (1 -γ d )] p d γ (1 -ξ d ) > 0,
where γ is defined in Proposition 2. Moreover, for any u ∈ P, we have whp

lim inf C(u, q) ñ ≥ lim | P(n) | ñ > 0 (ii) If q > q c
, for an uniformly chosen player u, we have C(u, q) = o p (ñ). The same result holds if o(n) players are chosen uniformly at random.

When γ r = 0 for all r ≥ 0, we recover the result of [START_REF] Lelarge | Diffusion and cascading behavior in random networks[END_REF]. If γ r = 1 for all r ≥ 1 (which means that we systematically replace each vertex i of degree d i by a clique of size d i ) we can even be more precise in the case q < q c : we have in this case, C(u, q) = | P(n) | for any u ∈ P(n) since in this case, the contagion starting from a pivotal player will propagate to the pivotal players only.

From these results, we see that the impact of clustering is different for low values of the mean degree and for high values of the mean degree. In the low values regime, as the clustering increases, the contagion threshold decreases whereas in the high values regime, the opposite happens. We see that in the low values regime for the mean degree, the clustering makes the contagion more difficult whereas it 'helps' the contagion in the high values regime.

IV. DIFFUSION THRESHOLD FOR RANDOM GRAPHS WITH CLUSTERING

A. Diffusion model

In this section, we study a simple diffusion model depending on a parameter π ∈ [0, 1] which can be described in term of a bond percolation process in a general graph G. Randomly delete each edge with probability 1π independently of all other edges. denote by G π the resulting graph. Then any active node will activate all nodes in its component in G π . As in previous section, we will derive conditions under which a single starting active node can activate a large fraction of the population in G = G (n, d, γ). This problem corresponds to the existence of a 'giant component' in the random graph obtained after bond percolation. Note that this model of diffusion corresponds to a simple epidemics with probability of contagion given by the parameter π ∈ [0, 1].

B. Phase transition for the diffusion

In order to state our result, we first need to recall some basic results about random graphs with small order. For d ∈ N, let K d be the complete graph on d vertices denoted {1, . . . , d}, with d(d -1)/2 edges. For π ∈ [0, 1], we denote by K d (π) the random graph obtained from K d after bond percolation with parameter π, i.e. each edge of K d is kept independently of the others with probability π, otherwise it is removed.

We need to compute the probability that the component in K d (π) containing vertex 1 has k vertices denoted by f (d, k, π). Note that f (d, d, π) is simply the probability that K d (π) is connected and has been computed in [START_REF] Gilbert | Random graphs[END_REF]. Indeed simple computations show that we have the simple recurrence relation

f (d, d, π) = 1 - d-1 k=1 d -1 k -1 f (k, k, π)(1 -π) k(d-k) , f (d, k, π) = d -1 k -1 f (k, k, π)(1 -π) k(d-k) , (1)
for any k ≤ d.

We now define for d ∈ N and π ∈ [0, 1], the random variable K(d, π, γ) defined by

P (K(d, π, γ) = k) = (1 -γ d )1(k = d) + γ d f (d, k, π),
where f is defined in [START_REF] Bollobás | Random graphs[END_REF].

For a graph G = (V, E) and a parameter ∈ [0, 1], we denote for any v ∈ V by C(v, π) the size of the component in the bond percolated graph G π containing v. . Let D * be a random variable with distribution p * r given by p * r-1 = rpr λ for all r ≥ 1. We define π c as the solution of the equation: Note that in the particular case where γ r = 0 for all r, we have K(d, π, 0) = d so that we get π c =

πE [K(D * + 1, π, γ) -1] = 1. • if π > π c , we have whp that lim inf C(u,π) ñ > 0, i.e

E[D]

E[D(D-1)] where D is the typical degree in the random graph and we recover a standard result in the random graphs literature (see Theorem 3.9 in [START_REF] Janson | On percolation in random graphs with given vertex degrees[END_REF]).

We can guess the value of the diffusion threshold π c using a branching process approximation. Indeed the random graph G (n, d) can be approximated by a branching process in which each node (except the root) has a number of offspring distributed as D * . The degree of a node v in the corresponding random tree is thus distributed as D * + 1. Let us assume D * + 1 = d. If we replace v by a clique K of size d with probability γ d , and delete independently each edge inside the clique with probability 1-π, then the probability that the component of v inside K contains k vertices is given by f (d, k, π). Hence the probability that v is linked to k vertices is:

(1 -γ d )1(d = k) + γ d f (d, k, π) = P (K(D * + 1, π, γ) = k).
The new distribution of offspring is thus K(D * + 1, π, γ) -1. Finally, we remove each edge with probability π, which gives πE [K(D * + 1, π, γ) -1] for the expected number of offspring.

For regular graphs, we obtain that the diffusion threshold increases as the clustering increases, as it was already observed in [START_REF] Gleeson | How clustering affects the bond percolation threshold in complex networks[END_REF].

V. FURTHER REMARKS AND CONCLUSIONS

If the characteristic parameter of the epidemic is above some threshold for the diffusion model, or under a certain threshold for the contagion model, then a global cascade is possible, starting from a single infected individual. An interesting question would be to study the effect of clustering on the cascade size, especially for the contagion model. The cascade size for the diffusion model is derived by heuristic means in [START_REF] Gleeson | Analytical results for bond percolation and k-core sizes on clustered networks[END_REF].

The diffusion threshold increases when the clustering increases (for random regular graphs), which makes the diffusion more difficult to spread. The effect of clustering on the contagion threshold depends on the value of the mean degree in the graph: for low values of the mean degree, we observe that clustering inhibits the contagion, and the contrary happens in the high values regime. One can wonder what is the effect of clustering on contagion threshold when degree distribution and degree-degree correlations are fixed: note that the simple consideration of random regular graphs does not provide significant information here, since increasing clustering in a random regular graph (adding cliques) does not change the contagion threshold with our graph model.

  the degree of a vertex chosen uniformly at random among the n vertices of G (n, d), and D a random variable with distribution p, (i) is equivalent to the fact that D n d -→ D (convergence in distribution).

Proposition 3 .

 3 where ñ is the number of vertices in the graph, including those of degree one or zero. The local clustering coefficient of vertices with degree one or zero is null by definition, hence the clustering coefficient in the graph can be very low if the graph contains a lot of such vertices, even if other vertices are highly clustered. To overcome this, we can consider another definition of clustering coefficient (see the next paragraph). a) Computation of the biased clustering coefficient: We consider the model G (n, d, γ) for a sequence d satisfying Condition 1 with probability distribution p = (p r ) ∞ r=0 . Then we have for the biased clustering coefficient of G (n, d, γ):

  there is no clustering and as γ increases, the biased clustering coefficient also increases to 1 -2-(p1+2p0) λ . b) Another definition of clustering coefficient: We keep the notations of the previous paragraph: d v is the degree of vertex v, and P v is the number of pairs of neighbors of v that share an edge together. Then we define the clustering coefficient C (n) 2

Proposition 4 .

 4 We consider the model G (n, d, γ) for a sequence d satisfying Condition 1 with probability distribution p = (p r ) ∞ r=0 . Then we have for the clustering coefficient of G (n, d, γ):

Fig. 1 .Fig. 2 .

 12 Fig. 1. Correlation between the clustering coefficient C 2 and the mean degree λ in G (n, d, γ) (when the initial graph G (n, d) has power law degree distribution pr ∝ r -τ e -r/50 , and varying parameters are τ and γ)

Fig. 3 .

 3 Fig. 3. Evolution of the clustering coefficient C 2 with respect to the clique probability γ, when the mean degree λ in G (n, d, γ) is fixed (for initial graph G (n, d) distributed as Erdős-Renyi)

Theorem 7 .

 7 Consider the random graph G = G (n, d, γ) for a sequence d satisfying Condition 1 with probability distribution p = (p r ) ∞ r=0

  and let B d be the number of vertices with degree d that are replaced by a clique. Then B d follows a Binomial distribution with parameters (n d , γ d ), where n d is the number of vertices with degree d in G (n, d). By Condition 1-(i), we have: n d /n → p d , so that the Law of Large Numbers implies: B d /n → p γ d p d (which is still true if p d = 0).

  . there exists a 'giant component' in the percolated graph G π .

• if π < π c , for an uniformly chosen player u, we have C(u, π) = o p (ñ). The same result holds if o(n) players are chosen uniformly at random, i.e. there is no 'giant component' in the percolated graph G π .


Agnès Lagnoux-Renaudie
email: lagnoux@cict.fr

Effective branching splitting method under cost contraint

Keywords: splitting method, simulation, cost function, Laplace transform, Galton-Watson, branching processes, iterated functions, rare event MSC : 65U05, 44A10, 60J80

This paper deals with the splitting method first introduced in rare event analysis. In this technique, the sample paths are split into R multiple copies at various stages during the simulation. Given the cost, the optimization of the algorithm suggests to sample a number of subtrials which may be non-integer and even unknown but estimated. To avoid this problem, we present in this paper three different approaches which provide precise estimates of the relative error between P(A) and its estimator.

Introduction

General settings

The study of rare events is an important area in the analysis and prediction of major risks as earthquakes, floods, air collision risks, etc. Studying the major risks can be taken up by two main approaches which are the statistical analysis of collected data and the modelling of the processes leading to the accident. The statistical analysis of extreme values needs a long observation time since the very low probability of the events considered. The modelling approach consists first in formalizing the system considered and then in using mathematical (Aldous [START_REF] Aldous | Probability Approximations via the Poisson Clumping Heuristic[END_REF] and Sadowsky [START_REF] Sadowsky | On Monte Carlo estimation of large deviations probabilities[END_REF]) or simulation tools to obtain some estimates.

Analytical and numerical approaches are useful, but may require many simplifying assumptions. On the other hand, Monte Carlo simulation is a practical alternative when the analysis calls for fewer simplifying assumptions. Nevertheless, obtaining accurate estimates of rare event probabilities, say about 10 -9 to 10 -12 , using traditional techniques require a huge amount of computing time.

Many techniques for reducing the number of trials in Monte Carlo simulation have been proposed, the first one is based on importance sampling (IS), e.g. de Boer [START_REF] De Boer | Analysis and efficient simulation of queueing models of telecommunication systems[END_REF] and Heidelberger [START_REF] Heidelberger | Fast Simulation of Rare Events in Queueing and Reliability Models[END_REF]. Fundamentally, IS is based on the notion of modifying the underlying probability distribution in such a way that the rare event occurs much more frequently. But the use of IS requires a deep knowledge of the studied system. An alternative way is to use trajectory splitting, based on the idea that there exists some well identifiable intermediate system states that are visited much more often than the target states themselves and behave as gateway states to reach the rare event. In contrast to IS type algorithms, the step-by-step evolution of the system follows the original probability measure. Thus we consider a decreasing sequence of events B i leading to the rare event A :

A = B M +1 ⊂ B M ⊂ . . . ⊂ B 1 .

Then P(A) is given by P(A) = P(A|B M )P(B M |B M -1 ) . . . P(B 2 |B 1 )P(B 1 ), [START_REF] Aldous | Probability Approximations via the Poisson Clumping Heuristic[END_REF] where on the right hand side, each conditioning event is ≪ not rare ≫ . For the applications we have in mind, these conditionnal probabilities are in general not available explicitly. Instead we know how to make evolve the particles from level B i to the next level B i+1 (e.g. Markovian behavior).

The principle of the algorithm is at first to run simultaneously several particles starting from the level B i ; after a while, some of them have evolved ≪ badly ≫ , the other have ≪ well ≫ evolved i.e. have succeeded in reaching the threshold B i+1 . ≪ Bad ≫ particles are then moved to the position of the ≪ good ≫ ones and so on until A is reached. In such a way, the more promising particles are favoured ; unfortunately that algorithm is hard to analyse directly because of the interaction introduced between particles. Examples of this class of algorithms can be found in [START_REF] Aldous | ≪Go with the winners≫ algorithms[END_REF] with the ≪ go with the winners ≫ scheme, in [START_REF] Diaconis | Three examples of Monte-Carlo Markov chains : at the interface between statistical computing, computer science, and statistical mechanics, Discrete probability and algorithms[END_REF][START_REF] Jerrum | The Markov chain Monte Carlo method : an approach to approximate counting and integration[END_REF] in the context of the approximate counting and in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications, Probability and its Applications[END_REF][START_REF] Doucet | An introduction to sequential Monte Carlo methods, Sequential Monte Carlo methods in practice[END_REF] in a more general setting.

Nevertheless, in practice the trajectory splitting method may be difficult to apply. For example, the case of the estimation of the probability of a rare event in random dynamical systems is more complex, since the difficulty to find theoretically the optimal B i . Furthermore, the probability to reach B i varies generally with the state of entrance in level B i-1 . But it is not always the case e.g. for Markovian models (like diffusions).

A mathematical tool well adapted to study this type of algorithms is the Feynman-Kac approach developped in [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications, Probability and its Applications[END_REF]. Asymptotic results are derived, such as LLN, CLT, and Large Deviations principles ; in particular asymptotic variance of the estimator of the rare event probability is given. Non asymptotic results such as uniform Lp mean error bounds and exponential concentration inequalities with respect to the time horizon can be also found in this relevant book. Getting precise confidence intervals is more challenging. Nevertheless, all these algorithms lie on a common base, simpler to analyse and called branching splitting model. In this technique, interactions between particles are avoided and its relative simplicity allow us to derive precise results (Chernoff type bound of the estimator) and to have better knowledge and understanding on splitting models in general. We must precise here that we consider only one dimensional models as introduced in Garvels [START_REF] Garvels | The splitting method in rare event simulation[END_REF] or in a more refined version : the RESTART method [START_REF] Villén-Altamirano | RESTART : a Method for Accelerating Rare Event Simulations[END_REF][START_REF] Villén-Altamirano | Restart : An Efficient and General Method for Fast Simulation of Rare Event[END_REF].

In the branching splitting technique, make a {0, 1} Bernoulli trial to check whether or not the set event B 1 has occured. In that case, we split this trial in R 1 Bernoulli subtrials, and for each of them we check again whether or not the event B 2 has occured. This procedure is repeated at each level, until A is reached. If an event level is not reached, neither is A , then we stop the current retrial. Using N independent replications of this procedure, we have then considered NR 1 . . . R M trials, taking into account for example, that if we have failed to reach a level B i at the i-th step, the R i . . . R M possible retrials have failed. Clearly the particles reproduce and evolve independently.

An unbiased estimator of P(A) is given by the quantity

P = N A N M i=1 R i , (2)
where N A is the total number of trajectories having reached the set A. Considering that this algorithm is represented by N independant Galton-Watson branching processes (Z n ) n , as done in [START_REF] Lagnoux | Rare event simulation[END_REF], the variance of P can be then derived and depends on the probability transitions and on the mean numbers (m i ) of particles successes at each level. Lead by the heuristic presented in [START_REF] Villén-Altamirano | RESTART : a Method for Accelerating Rare Event Simulations[END_REF][START_REF] Villén-Altamirano | Restart : An Efficient and General Method for Fast Simulation of Rare Event[END_REF], an optimal algorithm is derived by minimizing the variance of the estimator for a given budget (computational cost), defined as the expected number of trials generated during the simulation, where each trial is weighted by a cost function.

The optimization of the algorithm (Lagnoux [START_REF] Lagnoux | Rare event simulation[END_REF]) suggests to take all the transition probability equal to a constant P 0 and the number of splitting equal to the inverse of this constant. We then deduce the number of thresholds M and finally N is given by the cost. This result is not surprising since it means that the branching processes are critical Galton-Watson processes. In other words, optimal values are chosen in such a way to balance the loss of variance from too little splitting and the exponential growth in computational effort from too much splitting. Now, we are interested in the study of the precision of this algorithm by deriving an upper bound to the quantity

P(| P -P(A)|/P(A) ≥ α) (3)
The Chebycheff's bound being too crude, we will use the Chernoff type bound based on the Laplace transfom of the normalized Galton-Watson branching process

W M +1 := Z M +1 /E(Z M +1 ).
We therefore need to get estimates on the Laplace transform of W M +1 which depends on the n-th iterate of a function ψ. In practice, the optimal number of thresholds M is not very large and therefore asymptotic estimates may not be accurate. That optimal number being derived from these bounds, a numerical approach would be unworth too. Hence we want to derive explicit upper and lower bounds for a given number of thresholds. For an asymptotic analysis of the branching splitting model, the reader is referred to Glasserman et al. [START_REF] Glasserman | A large deviations perspective on the efficiency of multilevel splitting[END_REF][START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF].

In this paper, precise estimates are derived using the following technique : instead of a single bounding function for ψ, several ones are used to obtain sharper upper and lower estimates. These bounding functions are chosen in the low dimensional Lie groups of the homographic and affine functions. The higher the dimension of the Lie group, the more precise is the approximation since the dimension describes the number of parameters to adjust the function with. Unfortunately, we did not find higher dimensional Lie algebras of monotone functions (necessary property to iterate merely the inequalities obtained). The interest of using such functions lies also in the fact that their iterates can be explicitly computed. This technique leads to accurate bounds on the probability (3) (see Proposition 1).

In practice, the adjustment of P i close to the optimal value may be done during a first phase. The proportion of the cost devoted to this learning part will be the topic of a forthcoming paper (see Section 7). But it soon appears that, even in the case of P i 's close to optimals, the fact that the number of replicas is not an integer destroys rapidly the accuracy of the algorithm : in such a case, one can take R i equal to the closest integer (k or k + 1) of the optimal value R but whatever the choice we have made, the criticality of the Galton-Watson process will be lost and the loss of precision is significant, see Figure 1 and Tabular 1. Remark that Random 1 provides the closest results from the optimals.

P(A) = 3.
R i = [R] + 1 R-[R]≥0.

The results

Using the technique based on low dimensional Lie groups, described in section 1.1, we shall prove the following precise estimates :

Proposition 1 For a given cost C, there exists a generic constant α 1 such that for α ≤ α 1 ,

P | P -P(A)| P(A) ≥ α ≤      2 exp -Cα 2 P 0 8(1-P 0 )h(P 0 ) m M m-1 m M +1 -1 2 for m = 1 2 exp -Cα 2 P 0 8(1-P 0 )h(P 0 ) 1 (M +1) 2 for m = 1
The value of α 1 will be given in Section 3.

Remark 2 For large value of M, the bound of (3) then behaves like 2 exp -

Cα 2 P 0 8(1 -P 0 )h(P 0 ) 1 m M +2 for m > 1 2 exp - Cα 2 P 0 8(1 -P 0 )h(P 0 ) 1 (M + 1) 2 for m = 1 2 exp - Cα 2 P 0 8(1 -P 0 )h(P 0 ) m M for m < 1
With ( 13) and [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF] in Section 2, we obtain the following gaussian bound

2 exp - α 2 8 N var(W M +1 ) (4)
Considering R as a random variable in the two ways described in Section 1.1 and using the same type of techniques we get Proposition 3 (1) Sampling at each success There exists generic constants c 2 and α 2 such that for α ≤ α 2

P | P -P(A)| P(A) ≥ α ≤ 2 exp -c 2 C α 2 h(P 0 ) 1 (M + 1) 2
where c 2 is a constant depending on p and k.

(

2) Sampling a random environment

There exists a generic constant α 3 such that for α ≤ α 3

P | P -P(A)| P(A) ≥ α ≤ 2E exp - α 2 8(1 -P 0 ) 1 M i=0 ξ -1 0 . . . ξ -1 i C where ξ 0 = P 0 , ξ i = R i P 0 for i = 1 . . . M, C = C h(P 0 ) E(ξ)-1 E(ξ) M +1 -1 .
The paper is organized as follows. Section 2 recall quickly the general settings of the branching splitting model. We obtain, in Section 3, estimates on (3) and their sensitivity in m. In Sections 4 and 5, the two random models are studied and we derive precise estimates on (3). Section 6 presents a numerical illustration. Finally in Section 7, we conclude and discuss the merits of this approach and potential directions for further researches.

As said in Introduction and following [START_REF] Lagnoux | Rare event simulation[END_REF], we consider N independent Galton-Watson branching processes (Z (i) n ) n≥0 , i = 1, . . . , N where for each i, Z (i) n is the number of particles derived from the i-th particle (Z (i) 0 =1) that have reached the level B n . Then, letting R i the sampling number at level i,

P := 1 N N i=1 P i , where P i = Z (i) M +1 R 1 . . . R M . (5)
To lighten notation, we will consider only the case N = 1 in the following, i.e. we will consider the process (Z n ) n≥0 with Z 0 = 1. We have the following recurrence relation

Z n+1 = Zn j=1 X (j) n (6)
where for each n, the random variables (X (j) n ) j≥1 are i.i.d. with common law a Binomial distribution with parameters (R n , P n+1 ) for n ≥ 1 and a Bernoulli distribution with parameter P 1 for n = 0. Let introduce the following quantities

m 0 = P 1 , m n = R n P n+1 , n = 1, . . . , M + 1,
which are the mean number of particles success at each level. Then

P = Z M +1 R 1 . . . R M = P(A) Z M +1 m 0 . . . m M (7)
According to [START_REF] Babillot | The random difference equation X n = A n X n-1 + B n in the critical case[END_REF], we obtain that E(Z n+1 |Z n ) = m n Z n , so using this relation repeatedly, we find that E(Z n+1 ) = n i=0 m i . Introducing the random variable W n = Z n /E(Z n ) gives the new expression

P = P(A)W M +1 (8)
from which, we easily deduce that the estimate P is unbiased. Integrating N, the variance of P is given by var

( P ) = P(A) 2 N M i=0 1 -P i+1 m 0 . . . m i . ( 9
)
The minimization of the variance of P(A) for a given budget C defined by

C = N h(P 1 ) + M i=1 h(P i+1 )m 0 . . . m i-1 R i (10)
leads to the optimal parameters of the algorithm given by [START_REF] Lagnoux | Rare event simulation[END_REF]

             P i = P 0 , i = 1, . . . , M + 1, m i = 1, i = 1, . . . , M, N = C/[(M + 1)h(P 0 )], M = ⌊log P(A)/y 0 ⌋,
where y 0 is the solution of some equation which depends on the unit cost function h. It means that the branching processes are critical Galton Watson processes. The reader is referred to Harris [START_REF] Harris | The theory of branching processes, Dover Phoenix Editions[END_REF], Lyons [START_REF] Lyons | Probability on trees and networks[END_REF] and Athreya and Ney [START_REF] Athreya | Branching processes, Die Grundlehren der mathematischen Wissenschaften[END_REF] for more details on Galton-Watson and branching processes.

Remind the goal of this paper is to study the precision of the algorithm by deriving an upper bound to (3) that decomposes in two probabilities itselves bounded by

P P ≥ (1 + α)P(A) ≤ exp N inf u>0 log E e uW M +1 -u(1 + α) , (11)
P P ≤ (1 -α)P(A) ≤ exp N inf u<0 log E e uW M +1 -u(1 -α) . ( 12
)
Thus precise estimates of the Laplace transform of W M +1 need to be derived.

In the following, we consider our model with all the P i 's equal to a constant P 0 and all the R i 's equal to R ; but m = RP 0 could be different from the optimal value 1. From ( 9) and ( 10), remark

var(W M +1 ) := σ 2 M +1 =    1 P 0 -1 m M +1 -1 (m-1)m M for m = 1, 1 P 0 -1 (M + 1) for m = 1, (13)
and

C =    Nh(P 0 ) m M +1 -1 m-1 for m = 1, Nh(P 0 )(M + 1) for m = 1. ( 14
)
3 First model : taking R deterministic with m = 1

The model

Throughout this section, we consider the branching process associated to the splitting algorithm in the case R deterministic with m = 1 and we aim at deriving relevant confidence intervals.

Note the generating function of Z M +1 is P 0 f oM (s) + 1 -P 0 where f (s) = (P 0 s + 1 -P 0 ) R . To lighten notation, we denote the n-th functional iterate of a function f by f n .

Letting ψ (u) := log (f (e u )) = R log (P 0 e u + 1 -P 0 ) and h(u) = exp(u), we observe that ψ = h -1 of oh ; thus ψ M +1 = h -1 of M +1 oh and the Laplace transform of W M +1 becomes

E e uW M +1 = P 0 f M e u P 0 m M + 1 -P 0 = exp 1 R ψ M +1 u P 0 m M (15)
As a consequence, ( 11) and ( 12) yield to

P | P -P(A)| P(A) ≥ α ≤ exp N inf u>0 F + M +1 (u) + exp N inf u<0 F - M +1 (u)
where

   F + M +1 (u) := 1 R ψ M +1 u P 0 m M -(1 + α) u F - M +1 (u) := 1 R ψ M +1 u P 0 m M -(1 -α) u
Remark 4 Note that for u ≤ 0, E e uWn ≤ 1 and so E e uW ≤ 1. What happens for u ≥ 0 ? We have

ψ n+1 u P 0 m n = R log E e uW n+1
As a consequence, information on ψ n+1 u P 0 m n provides us results on the exponential integrability of W n+1 and W .

The Laplace transform of W n+1

Introduction

Wanting to obtain precise bounds on the Laplace transform of W n+1 given by (15) (M being replaced by n), we are interested in the behavior of ψ n+1 . But since ψ n+1 has no explicit expression, we are brought to estimate ψ in a first step and then to iterate these estimates.

Heuristic for estimating ψ : instead of using a single function to bound ψ, we use several bounding functions to obtain sharper upper and lower estimates. Near the origin, we will use homographic functions and elsewhere affine functions.

For example, for u ≥ 0 and m > 1, near the origin a ≪ good ≫ homographic function (same value at 0, same first and second derivatives at 0) gives a sharp upper bound. For large u, we simply use the asymptotic direction. We get

ψ(u) ≤    mu/ 1 -1-P 0 2 u for 0 ≤ u ≤ u 0 Ru for u 0 ≤ u for u 0 such that mu 0 / 1 -1-P 0 2 u 0 = Ru 0 .
Let us introduce h(u; x, y) = xu 1-yu and g(u; x, y) = xu + y. More generally, we obtain upper and lower estimates of the following form

   h(u; x, y) for u close to 0 g(u; x, y) otherwise
We choose these bounding functions in the low dimensional Lie groups of the homographic and affine functions. The higher the dimension of the Lie group, the more precise is the approximation since the dimension describes the number of parameters to adjust the function with (for example, the k first derivatives at a given point). Unfortunately, we did not found higher dimensional Lie algebras of functions having the monotonicity assumption required to iterate merely the inequalities obtained for ψ. The interest of using such functions lies also in the fact that their iterates can be explicitly computed. Obviously

h n (u; x, y) = x n u 1 -yu x n -1 x-1 and g n (u; x, y) = x n u + y x n -1 x -1
The problem is that the set formed by the homographic functions and the affine functions is not a group ; so a crucial point will be to determine the number of iterations corresponding to a change of group.

Heuristic for estimating ψ n : Since the estimates of ψ have two expressions, we will obtain estimates of ψ n with three functions. Indeed, starting from u close to 0, we iterate h, starting from u far from 0, we iterate the affine function g. Nevertheless when u starts from a midle value, we shall encounter a change of regime as described below.

More precisely, in the case when m > 1 and u > 0 e.g., 0 is a repulsive fix point (see next section). So starting from u ≥ 0, ψ n (u) diverges while n goes to ∞ and there are different regimes for the upper estimates :

-if u ≥ u 0 , we simply iterate u → Ru and the n-fold convolution of ψ is lower than the n-fold convolution of Ru, -if u ≤ u 0 , we start by iterating h .; m, 1-P 0 2 and if fortunately after n iterations h n u; m, 1-P 0 2 is still lower than u 0 then the n-fold convolution of ψ is lower than the n-fold convolution of h. In the other case, after k iterations, h k u; m, 1-P 0 2 is greater than u 0 then we continue iterating u → Ru. The value k 0 appearing in the following lemma corresponds to that change of regime.

Finally, the estimates of ψ n reflects that change of regime and have the following form

      
h n (u; x, y) for u close to 0 g n (u; x, y) for u far from 0 au α otherwise

Behavior and estimates of ψ for m > 1

We plot in Figure 2 ψ in the case m > 1. Note that ψ has two fix points : 0 which is repulsive and some q < 0 which is attractive (s.t. ψ(q) = q). Proceeding as explained in the previous section, straightforward studies of Behavior of ψ for m>1 0 is a repulsive fix point q is an attractive fix point

ψ (u) ≤    h u; m, 1-P 0 2 for 0 ≤ u ≤ u 0 Ru for u ≥ u 0 ψ (u) ≥    mu for u ≤ v 0 g (u; R, R log (P 0 )) for u ≥ v 0
where u 0 := 2 and v 0 :=log (P 0 ) /(1 -P 0 ) ≥ 0.

Lemma 3.2 Estimates of ψ for u ≤ 0 and m > 1. Let -b such that h (b; m, 1 -P 0 ) < b, -d := [q -h (b; m, 1 -P 0 )]/(q -b), -m q := ψ ′ (q) = me q( 1-1 R ) , -u -∞ := R log (1 -P 0 ), -v 1 := arg u {h (u; m, 1 -P 0 ) = g (u; m q, q(1 -m q))}, -v 2 := arg u {g (u; m q, q(1 -m q)) = u -∞ } = q + (u -∞ -q)/m q. ψ (u) ≤          h (b; m, 1 -P 0 ) for b ≤ u ≤ 0 g (u; d, q(1 -d)) for q ≤ u ≤ b q + h u -q; m q, (1 -P 0 )e -q R for u 2 ≤ u ≤ q for some u 2 ≤ 0 ψ (u) ≥        h(u; m, 1-P 0 2 ) for v 1 ≤ u ≤ 0 g (u; m q, q(1 -m q)) for v 2 ≤ u ≤ v 1 u -∞ for u ≤ v 2 3.2.3 Estimates of ψ n for m > 1 Theorem 3.1 Estimates of ψ n for u ≥ 0 and m > 1 Let u ⋆,n := h -1 n u 0 ; m, 1-P 0 2 = 2/[m n (1 + σ 2 n /R)] and γ := log R/ log m > 1.
Then there exists δ i for i = 0, 1 such that

ψ n (u) ≤              h n (u; m, 1-P 0 2 ) for 0 ≤ u ≤ u ⋆,n u 0 R n-δ 0 u 1+ 1-P 0 2(m-1) u γ for u ⋆,n ≤ u ≤ u 0 R n u for u 0 ≤ u ψ n (u) ≥        m n u for 0 ≤ u ≤ v 0 m n R n-δ 1 u γ (v 0 -q g ) for v 0 m n ≤ u ≤ v 0 g n (u; R, R log P 0 ) for v 0 ≤ u Proof Upper bound During the proof we write h instead of h(.; m, 1-P 0 2 ). If u ≤ u ⋆,n , -since m > 1, u ≤ h (u) and since h is increasing, u ≤ h (u) ≤ . . . ≤ h n (u). -since u ≤ u ⋆,n , h n (u) ≤ u 0 then ψ n (u) = ψ n-1 (ψ (u)) ≤ ψ n-1 (h (u)) ≤ . . . ≤ h n (u) = m n u 1 -1-P 0 2 u m n -1 m-1 If u ≥ u 0 , since R ≥ 1, u 0 ≤ u ≤ Ru ≤ . . . ≤ R n u and so ψ n (u) ≤ R n u If u ∈ [u ⋆,n , u 0 ]
, since u ≤ u 0 , we start by iterating h and since u ≤ h (u) and u ≥ h -1 n (u 0 ), the iterates h k (u) finally become greater than u 0 and then we iterate u → Ru. Let k 0 such that h k 0 (u) ≤ u 0 and h k 0 +1 (u) > u 0 . k 0 corresponds to the change of regime and

k 0 := log   u 0 1 + 1-P 0 2(m-1) u 0   / log m -log   u 1 + 1-P 0 2(m-1) u   / log m Letting δ 0 = log u 0 1+ 1-P 0 2(m-1) u 0 / log m, we get ψ n (u) = ψ n-k 0 (ψ k 0 (u)) ≤ R n-k 0 h k 0 (u) ≤ R n-k 0 u 0 ≤ . . . ≤ u 0 R n-δ 0   u 1 + 1-P 0 2(m-1) u   γ
A similar behavior is observed in the other cases thus we will omit their proofs.

Theorem 3.1 allows us to state about the exponential integrabibilty of W , the limit of W n as n → +∞ :

Corollary 5 The random variable uW γ ′ is exponentially integrable, independently of the value of u, where

γ ′ = γ γ-1 ≥ 1 i.e. E e uW γ ′ < +∞ for all u ∈ R Theorem 3.2 Estimates of ψ n for u ≤ 0 and m > 1. Let -u ⋆⋆,n := h -1 n (b; m, 1 -P 0 ) = b/[m n (1 + σ 2 n b/R)], -v ⋆,n := h -1 n v 1 ; m, 1-P 0 2 = v 1 /[m n (1 + σ 2 n v 1 /(2R))], -v ⋆⋆,n := g -1 n (v 2 ; m q, q(1 -m q)) = q -(q -v 2 )/m n ,
η := log d/ log m and µ := log m/ log m q .

Then there exists δ i for i = 2, 3 such that

ψ n (u) ≤                  h n (u; m, 1 -P 0 ) for u ⋆⋆,n ≤ u ≤ 0 q + d n-δ 2 u 1+ 1-P 0 m-1 η (b -q) for b ≤ u ≤ u ⋆⋆,n g n (u; d, q(1 -d)) for q ≤ u ≤ b q + h n u -q; m q, (1 -P 0 )e -q R for u 2 ≤ u ≤ q ψ n (u) ≥                    h n (u; m, 1-P 0 2 ) for v ⋆,n ≤ u ≤ 0 q + m n-δ 3 q u 1+ 1-P 0 2(m-1) u -µ (v 1 -q) for v 1 ≤ v ⋆,n g n (u; m q, q(1 -m q )) for v ⋆⋆,n ≤ u ≤ v 1 u -∞ for u ≤ v ⋆⋆,n

Estimates of ψ n for m < 1

Proceeding in the same way, we obtain the same kind of inequalities.

Confidence intervals

≥ α

In both cases, for α small enough and m = 1,

P | P -P(A)| P(A) ≥ α ≤ h + + h - where        h + := exp - 2C h(P 0 )(1-P 0 )R m M +1 1-m 1-m M +1 2 √ 1 + α -1 2 h -:= exp - C h(P 0 )(1-P 0 )R m M +1 1-m 1-m M +1 2 1 - √ 1 -α 2
Proof All the cases behave similarly : so we just treat the case m > 1 and u ≥ 0 for example. We have

F + M +1 (u) := 1 R ψ M +1 u P 0 m M -(1 + α) u. Then inf u>0 F + M +1 (u) = inf inf u∈[0,u ⋆,M +1 ] F + M +1 (u) , inf u∈[u ⋆,M +1 ,u 0 ] F + M +1 (u) , inf u>u 0 F + M +1 (u) For u ∈ [0, u ⋆,M+1 ], F + M +1 (u) ≤ u 1 -σ 2 M +1 u/2 -u (1 + α) := φ M +1 (u)
The cancellation of the first derivative of φ M +1 gives us the minimum of φ M +1 :

u s = 2 σ 2 M +1 1 - 1 √ 1 + α
which is lower than u ⋆,M +1 for α small enough i.e.

α ≤ 1 -u ⋆,M +1 σ 2 M +1 /2 -2 -1 For u ∈ [u ⋆,M+1 , u 0 ], F + M +1 (u) ≤ R M -δ 0 u 0 u 1+ 1-P 0 2(m-1) u ⋆,M +1 P 0 m M γ -(1 + α) u := φ M +1 (u) . The cancellation of the first derivative of φ M +1 gives us the mini- mum of φ M +1 : ũs = (1 + α) γu 0 R M -δ 0 P 0 m M + 1 -P 0 2(m -1) u ⋆,M +1 γ 1 γ-1 which is lower than u ⋆,M +1 for α ≤ γu 0 R M -δ 0 u ⋆,M +1   u ⋆,M +1 P 0 m M + 1-P 0 2(m-1) u ⋆,M +1   γ - 1
And so the solution leaves the interval.

For u ≥ u 0 ,

F + M +1 (u) ≤ R M u/ P 0 m M -u (1 + α) = u 1/P M +1 0 -(1 + α) ,
and 1/P M +1 0 -(1 + α) is positive for M great enough ; so the minimum is atteined in u 0 and the solution leaves the interval. Finally,

P P ≥ P (1 + α) is bounded by exp N R F + M +1 (u s ) = exp - N R 1 1 -P 0 m M +1 1 -m 1 -m M +1 √ 1 + α -1
Theorem 3.4 For m = 1 and α small enough, we have

P | P -P(A)| P(A) ≥ α ≤ g + + g - where        g + := exp - 2C h(P 0 )(1-P 0 )R 1 (M +1) 2 √ 1 + α -1 2 g -:= exp - C h(P 0 )(1-P 0 )R 1 (M +1) 2 1 - √ 1 -α 2
We check that this bound is better than the one for m = 1.

The next corollary follows from Theorem 3.3 and provides us an exact upper bound of P | P -P(A)|/P(A) ≥ α involving the variance of W M +1 that we might compare to (4), the approximate one obtained in the Introduction.

Corollary 6 For α such as in Theorem 3.3,

P | P -P(A)| P(A) ≥ α ≤ 2 exp - α 2 4 N var(W M +1 ) 1 - α 2
Proof By Theorem 3.3, the definition of N and (13), we deduce that

P | P -P(A)|/P(A) ≥ α is bounded by exp - 2N var(W M +1 ) √ 1 + α -1 2 + exp - N var(W M +1 ) 1 - √ 1 -α 2 But for α ≤ 2 and since √ 1 + α ≤ 1 + α 2 and √ 1 -α ≤ 1 -α 2 , √ 1 + α -1 2 = α 2 2 + α + 2 √ 1 + α ≥ α 2 4 1 1 + α 2 ≥ α 2 4 1 - α 2 1 - √ 1 -α 2 = α 2 2 -α + 2 √ 1 -α ≥ α 2 4 1 + α 2 ≥ α 2 4 1 - α 2
that leads to the corollary.

Proposition 1 and Remark 2 now follow from the previous theorems, using the same kind of argument as in Corollary 6. We plot in Figure 3 the different bounds using :

Numerical illustration

-the variance, -the Laplace transform in the optimal case : m = 1, -the Laplace transform in the other case : m = 1, in two different cases -for P(A) = 5 10 -9 and C = 6 10 7 , then m=1.055 (supercritical case), -for P(A) = 10 -11 and C = 2 10 8 , then m=0.924 (subcritical case).

4 Second model : sampling at each success

The model

Let k ∈ N such that k = ⌊1/P 0 ⌋. We decide here to allow R to vary. More precisely, each time a particle reach a higher level, we generate a realization of a Bernoulli random variable R whose parameter p = 1 -q will be derived by the same optimization scheme as used previously. Let ξ = RP 0 .

Here the estimator chosen is given by :

P = P(A) N N i=1 W i M +1 = 1 N N i=1 Z i M +1 E(Z M +1 ) = P(A) N N i=1 Z i M +1 P 0 E(ξ) M

The variance and its optimization 4.2.1 Study of the variance

Proposition 7 P is trivially an unbiased estimator and its variance is given by

var( P ) = P(A) 2 N 1 P 0 -1 M i=0 1 E(ξ) i + var(ξ) P 0 E(ξ) M i=1 1 E(ξ) i ( 16
)
Proof First of all, var( P

) = P(A) 2 N 1 [P 0 E(ξ) M ] 2 var(Z M +1
). The calculation of var( P ) then amounts to the calculation of var(Z M +1 ).

Note that for all X random variable and F filtration, var(X) = var(E(X|F )) + E( var(X|F )) [START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF] Applying ( 17) to X = Z M +1 and F = σ(Z M ) and then to X = X M and F = σ(R M ), and since

X i k ∼ Bin(R k , P 0 ), we get var(Z M +1 ) = var( Z M i=1 X i M ) = var(Z M E(X M |Z M )) + E(Z M var(X M |Z M )) = E(X M ) 2 var(Z M ) + var(X M )E(Z M ) = E(E(X M |R M )) 2 var(Z M ) + var(E(X M |R M )E(Z M ) + E( var(X M |R M ))E(Z M ) = E(ξ) 2 var(Z M ) + var(ξ)E(Z M ) + (1 -P 0 )E(ξ)E(Z M ) By a recurrent descent, var(Z M +1 ) = E(ξ) 2M var(Z 1 )+ var(ξ) M k=1 E(Z k )E(ξ) 2(M -k) +(1-P 0 )E(ξ) M k=1 E(Z k )E(ξ) 2(M -k)
But var(Z 1 ) = P 0 (1 -P 0 ) and E(Z k ) = P 0 E(ξ) k-1 , and so we finally get the result.

Optimization of the parameters

As done in Lagnoux [START_REF] Lagnoux | Rare event simulation[END_REF], an optimal algorithm is chosen via the minimization of the variance of P for a given budget C, keeping the optimal values for M and P 0 . The (average) cost is now

C = Nh(P 0 ) M i=0 E(ξ) i + NP 0 h(1) M -1 i=0 E(ξ) i (18)
Neglecting the cost introduced by the generation of the random splitting numbers at each success, we assume, in the following,

C = Nh(P 0 ) M i=0 E(ξ) i , (h(1) ≪ h(P 0 )). ( 19
)
Suppose E(ξ) = 1, then the variance can be rewritten as

var( P ) = P(A) 2 N 1 P 0 -1 + var(ξ) P 0 E(ξ) M i=0 1 E(ξ) i - var(ξ) P 0 E(ξ) = P(A) 2 N 1 P 0 -1 + var(ξ) P 0 E(ξ) 1 E(ξ) M E(ξ) M +1 -1 E(ξ) -1 - var(ξ) P 0 E(ξ)
and the cost has the following form

C = Nh(P 0 ) E(ξ) M +1 -1 E(ξ) -1
The optimal value of N is given by the cost

N = C h(P 0 ) E(ξ) -1 E(ξ) M +1 -1
and we have to minimize the expression

var( P ) = P(A) 2 CP 0 h(P 0 ) E(ξ) M +1 -1 E(ξ) -1 1 -P 0 + var(ξ) E(ξ) 1 E(ξ) M E(ξ) M +1 -1 E(ξ) -1 - var(ξ) E(ξ)
whose principal term is

P(A) 2 CP 0 h(P 0 ) 1 E(ξ) M E(ξ) M +1 -1 E(ξ) -1 2 which is minimal for E(ξ) = 1.

Confidence interval

From now on, we take

m = f ′ (1) = E(ξ) = P 0 (k + q) = 1 which is equivalent to p = 1 -q = 1 -δ.
Here the generating function of the reproduction law is given by f (s) = (P 0 e u + 1 -P 0 ) k (P 0 qe u + 1 -P 0 q) := exp {ψ(log s)} Remark 8 That form reflects the fact that R is randomly chosen between k and k + 1 following the parameter q.

Since log(P 0 s + 1 -P 0 ) ≤ 1 k+q log f (s) and m = 1,

E e uW M +1 ≤ exp 1 k + q log f M +1 (e u P 0 ) = exp 1 k + q ψ M +1 u P 0
Using again the Chernoff's bounding method,

P | P -P(A)| P(A) ≥ α = exp N inf u>0 F + M +1 (u) + exp N inf u<0 F - M +1 (u)
where

   F + M +1 (u) := 1 k+q ψ M +1 u P 0 -(1 + α) u F - M +1 (u) := 1 k+q ψ M +1 u P 0 -(1 -α) u
Once again we are interested in the Laplace of W M +1 which leads to study the behavior of the iterate of the function ψ.

The Laplace transform of W n+1

Let β := [k(1 -P 0 ) + q(1 -P 0 q)] /(k + q). We will proceed as explained in Section 3.2.1.

Estimates of ψ

Lemma 4.1 Estimates of ψ for u ≥ 0 There exists u 0 > 0 such that

ψ (u) ≤    h u; 1, β 2 for 0 ≤ u ≤ u 0 (k + 1)u for u ≥ u 0 Lemma 4.2 Estimates of ψ for u ≤ 0 There exists u 0 < 0 such that ψ (u) ≤ h(u; 1, β) for u ≤ u 0 4.3.2 Estimates of ψ n Proposition 9 Estimates of ψ n for u ≥ 0 Let u ⋆ := h -1 n (u 0 ; 1, β 2 ). ψ n (u) ≤          h n u; 1, β 2 for 0 ≤ u ≤ u ⋆ u 0 (k + 1) n-2 β ( 1 u -1 u 0 ) for u ⋆ ≤ u ≤ u 0 (k + 1) n u for u 0 ≤ u Proposition 10 Estimates of ψ n for u ≤ 0 ψ n (u) ≤ h n (u; 1, β) for u ≤ u 0

Confidence interval

Finally, as in the deterministic case, we deduce Theorem 4.1 For α small enough, we have

P | P -P(A)| P(A) ≥ α ≤ h + + h - (20)
where

       h + := exp - 2C (k+q)β[h(P 0 )+P 0 h(1)M/(M +1)] 1 (M +1) 2 √ 1 + α -1 2 h -:= exp - C (k+q)β[h(P 0 )+P 0 h(1)M/(M +1)] 1 (M +1) 2 1 - √ 1 -α 2
The first part of Proposition 3 follows from Theorem 4.1 using the same kind of argument as in Corollary 6 and provides better results than the results of the previous section since here we obtain an upper bound in exp{ 1 (M +1) 2 } like in the optimal case. 5 Third model : sampling a random environment

The model

Here we decide to sample a random environment described by (R 1 , R 2 , . . . , R M ) at the beginning of the simulation, to keep the optimal values of P i and M and to optimize the algorithm in E(R). More precisely, let r i = {R (i) j , j = 1 . . . M}, for i = 1 . . . N N sequence of i.i.d. Bernoulli random variables {k, k + 1}. r i is said to be the i-th ≪ environmental ≫ sequence. Given that environmental sequences, we consider N independent branching processes (Z (i) n ) n≥0 , i = 1, . . . , N in the same way as in the deterministic case except that each R j is a random variable associated to a random generating function (g.f.) f j (s) which represents the offspring's distributions in j-th generation with f j (s) = (P 0 s + 1 -P 0 ) R j .

(Z n ) is now a Branching Process in Randon Environment (BPRE).

As in usual case in random environments, we should be careful to distinguish the alea coming from the environment and the one from the process itself.

Since to each (Z (i) n ) n≤0 i = 1 . . . N of the N > 1 initial particles we associate a random environment r i , the N initial particles reproduces independently one from each other. As a consequence, we can rewrite P the estimator of P as the sum of N independent branching processes :

P = 1 N N i=1 Pi = P(A) N N i=1 Z (i) M +1 ξ 0 ξ 1 . . . ξ M = P(A) N N i=1 W (i) M +1
where ξ i = R i P 0 is a random variable for all i = 1 . . . M ; and we aim at minimizing its variance for a fixed effort to derive the optimal parameter for R.

Some general results on BPRE

For more general background and details on BPRE, see for instance Athreya and Karlin [START_REF] Athreya | On branching processes with random environments. I. Extinction probabilities[END_REF][START_REF] Athreya | Branching processes with random environments. II. Limit theorems[END_REF], Guivarc'h et al. [START_REF] Guivarc | Normalisation d'un processus de branchement critique dans un environnement aléatoire[END_REF], Smith and Wilkinson [START_REF] Smith | On branching processes in random environments[END_REF] and Tanny [START_REF] Tanny | Limit theorems for branching processes in a random environment[END_REF]. Suppose now that r = {r j , j ≥ 0} is a general sequence of i.i.d. random variables that determines the succession of offspring g.f.'s {f j (r j ; s), j ≥ 0},i.i.d., in a BPRE. Let us recall that, in complete analogy with the classical Galton-Watson process, a Galton-Watson process in random environment is sub-

critical if E (log f ′ 0 (1)) < 1, critical if E (log f ′ 0 (1)) = 1 and supercritical if E (log f ′ 0 (1)) > 1.
Moreover, in the subcritical and the critical (resp. supercritical) cases the probability of extinction given the environment r is one (resp. q(r) ∈ [0, 1[). Nevertheless in the critical case, under the assumption of a finite third moment for the random variables X n , the asymptotic conditional distribution of Z n /E (Z n |Z n > 0), given that Z n = 0, is exponential, expressing the extreme character of this event for large n. In the supercritical case, W n converges to some random variable W non-degenerated.

The variance and its optimization

Study of the variance

Proposition 11 P is trivially an unbiased estimator and its variance is given by var

( P ) = P(A) 2 N 1 P 0 -1 M i=0 E 1 ξ i ( 21
)
Proof First of all, var( P ) = P(A) 2 N var(W M +1 ). The calculation of var( P ) then amounts to the calculation of var(W M +1 ).

Applying [START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF]

to X = W M +1 = Z M +1 P 0 M i=1 ξ i and F = σ(Z M , r) and since W M is a martingale and X i M ∼ Bin(R M , P 0 ), we get var(W M +1 ) = var(W M ) + E 1 [P 0 M i=1 ξ i ] 2 var(Z M +1 |σ(Z M , r)) = var(W M ) + E Z M [P 0 M i=1 ξ i ] 2 var(X M |σ(Z M , r)) = var(W M ) + E Z M [P 0 M i=1 ξ i ] 2 R M P 0 (1 -P 0 ) = var(W M ) + (1 -P 0 )E W M P 0 M i=1 ξ i By a recurrent descent, var(W M +1 ) = var(W 0 ) + (1 -P 0 ) M k=0 E W k P 0 k i=1 ξ i
It remains to compute the expectation of

W k P 0 k i=1 ξ i
which is derived by induction and at step k we have

E W k P 0 k i=1 ξ i = 1 P 0 E 1 ξ k ( 22
)
Clearly the formula holds for k = 0 : E W 0 P 0 = 1 P 0 . To go from k -1 to k, assume [START_REF] Kesten | Random difference equations and renewal theory for products of random matrices[END_REF] for k -1.

E W

k P 0 k i=1 ξ i = E Z k-1 [P 0 k-1 i=1 ξ i ] 2 ξ k E(X k-1 |σ(Z k-1 , r)) = E W k-1 P 0 k i=1 ξ i = E W k-1 P 0 k-1 i=1 ξ i E 1 ξ
since the R i 's are i.i.d. And we finally get the result using the induction's hypothesis.

Optimization of the parameters

As done in Lagnoux [START_REF] Lagnoux | Rare event simulation[END_REF], an optimal algorithm is chosen via the minimization of the variance of P for a given budget C, keeping the optimal values for M and P 0 . The (average) cost is now

C = N h(P 0 ) M i=0 E(ξ) i + Mh(1) (23)
Neglecting the cost introduced by the generation of the random environment, we assume, in the following,

C = N h(P 0 ) M i=0 E(ξ) i , (h(1) ≪ h(P 0 )). (24)
Once the trivial cases have been isolated, we can suppose E(R) = 1 P 0 and E( 1 R ) = P 0 , then the variance can be rewritten as

var( P ) = P(A) 2 N ( 1 P 0 -1) E( 1 ξ ) M +1 -1 E( 1 ξ ) -1
and the cost has the following form

C = Nh(P 0 ) E(ξ) M +1 -1 E(ξ) -1
The optimal value of N is given by the cost

N = C h(P 0 ) E(ξ) -1 E(ξ) M +1 -1
and we have to minimize the expression

var( P ) = P(A) 2 C h(P 0 ) 1 P 0 -1 E( 1 ξ ) M +1 -1 E( 1 ξ ) -1 E(ξ) M +1 -1 E(ξ) -1 under the constraint E(R) = 2k +1-k(k +1)E( 1 R ).
To lead the analytic study, suppose that 1/P 0 = k + δ with δ ∈]0, 1[ and let

-u = E(ξ) -1 and v = E( 1 ξ ) -1, -ρ = (k+δ) 2 k(k+1) and α = δ(1-δ) k(k+1) .
The constraint becomes v = α -ρu and by Lagrange multipliers, we finally need to solve

F (M, u) = ρF (M, v) (25)
where

F (M, u) = (1+u) M (M u-1)+1 u[(1+u) M +1 -1] .
Remark 12 Note that u and v depend on M and lie in the following intervals

[-δ k+δ , 1-δ k+δ ] and [-1-δ k+δ , δ k ] respectively.
We may only state asymptotic results.

Proposition 13 Asymptotically in M, we only have three solutions :

-u → u 1 = 0 and v → v 1 = α -u → u 2 = α ρ and v → v 2 = 0 -u → u 3 = 1+α-ρ 2ρ and v → v 3 = α+ρ-1 2
Proof • Suppose first that u and v do not converge to 0 when M → +∞. So we have three cases to analyse (the case u < 0, v < 0 is not worth to consider since E(1/R)E(R) ≥ 1) : -u < 0 and v > 0 -u > 0 and v < 0 -u > 0 and v > 0 Case 1 : If u < 0 and v > 0 (or symetrically u > 0 and v < 0), since u does not converge to 0, there exists a subsequence of u such that u < -ǫ ∀M. Then asymptotically in M, (25) is equivalent to

- 1 u = ρ M 1 + v
which is absurd since v = -1 and u does not converge to 0. Case 2 : If u > 0 and v > 0, in the same way, there exists a subsequence of u and a subsequence of v such that u > ǫ and v > ǫ. Then asymptotically in M, (25) is equivalent to

M 1 + u = ρ M 1 + v
And given the constraint, u must necessarily converge to 1+α-ρ 2ρ and v to α+ρ-1 2 . Since the sequence u has only one adherence value and lies in a compact, it converges to u 3 = 1+α-ρ 2ρ and similarly v converges to v 3 = α+ρ-1 2 . To guarantee u > 0 and v > 0, δ must be in

[ k 1+2k , 1 2 ] := [δ 1 , δ 2 ].
• Suppose now that u converges to 0 when M → +∞. We have two cases to analyse : Mu bounded or not.

If Mu is not bounded, we can extract from Mu a subsequence which diverges. Then

F (M, u) = 1 -1 M u + 1 M u(1+u) M 1 - 1 (1+u) M +1 1 u(1 + u) ∼ 1 u(1 + u)
and thus ( 25) is equivalent to

1 u(1 + u) = ρ M 1 + v which is absurd.
So Mu is bounded and we can extract of Mu a subsequence which converges to 0 or to a constant C 1 . In the first case,

F (M, u) ∼ e M u (Mu -1) + 1 u(e M u -1) ∼ (1 + Mu)(Mu -1) + 1 u(1 + Mu -1) ∼ M
and ( 25) is equivalent to

M = ρ M 1 + v which is absurd since v > -1. So the subsequence of Mu converges to a constant C 1 that solves e C 1 (C 1 -1) + 1 C 1 (e C 1 -1) = ρ 1 + α Let f (x) = e x (x-1)+1
x(e x -1) -ρ 1+α . Note that f is strictly increasing, converges to 1 -ρ 1+α ≥ 0 when x → +∞ and to -ρ 1+α ≤ 0 when x → -∞ for δ ≤ δ 2 . So C 1 is defined uniquely and then the sequence Mu has only one limiting point and lies in a compact, so it converges to C 1 , u converges to u 1 = 0 and v converges to v 1 = α.

Symetrically we get the equivalent result for v : it converges to v 2 = 0 and u converges to u 2 = α ρ .

Remind that we want to determine the minimum of the variance of P .

Proposition 14 Let δ

0 := √ k( √ k + 1 - √ k). Then asymptotically in M, -if δ ∈ [0, δ 0 ], (u 1 , v 1 ) minimizes the variance, -if δ ∈ [δ 0 , 1], (u 2 , v 2 ) minimizes the variance. Proof Let β = P 2
C h(P 0 ) 1 P 0 -1 and for i = 1, 2, 3,

var i = β [1 + u i ] M +1 -1 u i [1 + v i ] M +1 -1 v i
First of all, compare var 1 and var 2 . For i = 1, 2,

var i ∼ β e (M +1)(1+C i -1/M ) -1 u i [1 + v i ] M +1 v i ∼ βM e C i -1 C i [1 + v i ] M +1 v i Thus var 1 var 2 ∼ e C 1 -1 C 1 C 2 e C 2 -1 u 2 v 1 1 + v 1 1 + u 2 M +1 But 1+v 1 1+u 2 < 1 ⇔ δ ≤ δ 0 := √ k( √ k + 1 - √ k). As a conclusion, var 1 ≤ var 2 ⇔ δ ≤ δ 0
Then note that var 3 is not defined on [0, δ 1 ] ∪ [δ 2 , 1] and δ 0 trivially belongs to ]δ 1 , δ 2 [. Consequently, for δ ∈ [δ 1 , δ 2 ], it is sufficient to compare -var 1 and var 3 on [δ 1 , δ 0 ], -var 2 and var 3 on [δ 0 , δ 2 ].

But

var 1 var 3 ∼ M e C -1 C u 3 v 3 α 4ρ 1+α (1+α+ρ) 2 ----→ M →∞ M +1 and 4ρ 1+α (1+α+ρ) 2 = 4 ρ 1+α (1+ ρ 1+α ) 2 < 1. Thus var 1 var 3 ----→ M →∞
0 and var 1 ≤ var 3 asymptotically in M. In the same way,

var 2 var 3 ----→ M →∞ 0.

Confidence interval

By analogy with the deterministic case, we are interested in bounds on

P | P -P(A)|/P(A) ≥ α So let ψ i (u) := R i log (P 0 e u + 1 -P 0 ), f 0,i := f 0 •. . .•f i and ψ 0,i := ψ 0 •. . .•ψ i ,
taking again ξ i = R i P 0 for all i = 1, . . . , M and ξ 0 = P 0 , and

ψ (2) (u) ≤    h u; m (2) , 1-P 0 2 for 0 ≤ u ≤ a + g(u; k + 1, b + -a + (k + 1)) for u ≥ a +
which leads to the following proposition Proposition 17 If a 0 = a + and for all i = 1 . . . M,

a i :=    a + if R i = k b + if R i = k + 1 ψ 0,M (u) ≤        h + 0,M (u) for 0 ≤ u ≤ u + g + M (u) for u + ≤ u ≤ a + M i=0 R i u + M i=0 i-1 j=0 R j a i -M i=0 i j=0 R j a + for u ≥ a + where -h + 0,M (u) = u M i=0 ξ i / 1 -1-P 0 2 u 1 + M i=1 ξ i . . . ξ M and -g + M (u) =    h + 0,M (u) w.p. P 0 ≤ h + 0,M (u) ≤ a + b + w.p. 1 -P (0 ≤ h 0,M (u) ≤ a + ) • For u ≤ 0. Let -a -:= arg u =0 {h(u; m (2) , 1 -P 0 ) = u} = 1-m (2)
1-P 0 , -b -:= h(a -; m (1) , 1 -P 0 ), -u -:= h -1 M +1 (a -; m (1) , 1 -P 0 ).

In the same way, we deduce ψ (1) (u) ≤    h u; m (1) , 1 -P 0 = m (1) u 1-(1-P 0 )u for a -≤ u ≤ 0 b -for u ≤ a - and ψ (2) (u) ≤    h u; m (2) , 1 -P 0 for a -≤ u ≤ 0 b -for u ≤ a - which leads to the following proposition Proposition 18

ψ 0,M (u) ≤        h - 0,M (u) for u -≤ u ≤ 0 g - M (u) for a -≤ u ≤ u - b -for u ≤ a -
where Subcritical case : E (log ξ) < 0 First of all, Diaconis and Freedman [START_REF] Diaconis | Iterated random functions[END_REF] allow us to state that the backward process ψ 0,M (u) = ψ 0 • . . . • ψ M (u) converges almost surely, at an exponential rate to a random limit that does not depend on the starting point u.

Concerning the random walk on the affine group in our particular case (ξ L = Ber on {kP 0 , (k + 1)P 0 }), we have by Maksimov [START_REF] Maksimov | A generalized Bernoulli scheme and its limit distributions[END_REF] y n ---→ But Kesten [START_REF] Kesten | Random difference equations and renewal theory for products of random matrices[END_REF] gives an asymptotic result on the distribution of D : for some C > 0 and some κ > 0, Critical case : E (log ξ) = 0 Maksimov [START_REF] Maksimov | A generalized Bernoulli scheme and its limit distributions[END_REF] asserts that the distribution of log y n / √ n approaches, as n → ∞, the truncated normal distribution for any starting point of the walk.

P(D ≥ x) ∼ x→∞ C x κ
Corollary 20 Let Z a truncated normal distribution and σ 2 = var(ξ).

P (0 ≤ h 0,n (u) ≤ a) ∼ n→∞ P Z ≤ - 1 √ n log u ∼ n→∞ 1 -e -log u 2 n 1 2
Supercritical case : E (log ξ) > 0 Following Maksimov [START_REF] Maksimov | A generalized Bernoulli scheme and its limit distributions[END_REF], the distribution of ( yn v ) 1/ √ n e -M √ n approaches, as n goes ∞, the log normal distribution for any starting point of the walk. in two different cases -for P(A) = 5 10 -9 and C = 6 10 7 , then m=1.055 (supercritical case), -for P(A) = 10 -11 and C = 2 10 8 , then m=0.924 (subcritical case). In this article, the relative simplicity of our model allows us to state explicit results as the Chernoff's bounds of the relative error between the estimate P and P(A). Going further in the calculus, using the heuristic presented above, one can also deduce a central limit theorem and so Berry-Esseen bounds. Then we study the sensitivity of the Chernoff's bounds depending on the choice of the splitting number R in three different algorithms based on the branching splitting model : when one can not be exactly in the critical case (which corresponds to the optimal algorithm), the best way to proceed is to consider R as a random variable that we generate at each success during the simulation. Besides, this procedure is currently used in practice, see for example [START_REF] Aldous | ≪Go with the winners≫ algorithms[END_REF].

In practice, we do not know the transition probablities but just empirical estimation on them, and we can bound to adjust the levels according to them. In a model where the P i are unknown but belong in some known interval, we may proceed in the following way :

(1) Choose an arbitrarly sequence R of (P i ) i=1...M +1 are derived.

Figure 1 .

 1 Figure 1. Sensitivity of the bounds in m for different values of α

  5 10 -11 and C = 10 8 P(A) = 5.10 -11 and C = 10 8 m = 1.111 Subcritical case : m = 0.909 Deterministic

  ......................................... ........... ........................................

Figure 2 .Lemma 3 . 1

 231 Figure 2. Behavior of ψ for m > 1

Theorem 3 . 3

 33 Bounds of P | P -P(A)| P(A)

Figure 3 .

 3 Figure 3. Upper estimates for m = 1

n→∞ v( 1

 1 + D)where D = ∞ i=0 ξ 0 . . . ξ i . In general, nothing is known on the distribution of D except that it satisfies the following functional equationF (λ) = pF λ kP 0 -1 (1 -p)F λ (k + 1)P 0 -1with p = P(R = k).

Corollary 19

 19 Letting κ =log(p(1 -p))/ log(k(k + 1)P 2 0 ), P (0 ≤ h 0,n (u) ≤ a) ∼ n→∞ P v(1 + D)

Figure 5

 5 Figure 5. Confidence interval

  . During the first step, sample a packet of θ 1,N N particles following the splitting algorithm with sampling numbers R (0) i . Thus empirical estimations P (1) i i=1...M +1

  = 3.5 10 -11 and C = 10 8 P(A) = 5.10 -11 and C = 10 8

		P(A) Variance	Laplace	Variance	Laplace
	Optimal m=1	0.05080	0.03948	0.04938	0.03835
		Supercritical case : m = 1.111 Subcritical case : m = 0.909
	Deterministic
	R i = [R] + 1 R-[R]≥0.5	0.05702	0.04428	0.05428	0.04215
	Random 1	0.05134	0.03990	0.04990	0.03878
	Random 2	0.05388	0.04026	0.05235	0.03910
	TAB2. Length of the 95% confidence interval given by the variance and the
		Laplace transform
	5	0.05702	0.04428	0.05428	0.04215
	TAB1. Length of the 95% confidence interval given by the variance and the
		Laplace transform
	This paper is concerned with the study of different ways to overcome this
	problem. Lead by [2], we choose at random the sampling number with the
	hope of improving the simulation. In a first model (Random1), we sample a
	Bernoulli random variable R i on {k, k+1} for each particle having reached level i started from level i-1 and we decide to adjust p := P(R 1 = k) such that m =
	1. A second model (Random2) consists in sampling a random environmental
	sequence (R 1 , R 2 , . . . , R M ) of M i.i.d. Bernoulli random variables R i on {k, k + 1} with common parameter p, derived by the same previous optimization
	approach with an additional constraint (the link between the expectations
	of R and its inverse). However, this problem is more complex and needs an
	approximate solution. Results are presented in Proposition 3 and in Tabular
	2.

and we get the result from[START_REF] Doucet | An introduction to sequential Monte Carlo methods, Sequential Monte Carlo methods in practice[END_REF].

Acknowledgement

I gratefully acknowledge my supervisors Dominique Bakry and Pascal Lezaud for their helpful suggestions and comments.

Finally,

And we are interested in the Laplace transform of W M +1 which leads to study the behavior of the iterate of the random functions ψ.

5.3

The Laplace transform of W n+1

Criticality

By simple arguments of convexity, Proposition 15 Asymptotically in M, -for all δ ∈ [0, δ 0 ], we are in the subcritical case, -for all δ ∈ [δ 0 , 1], we are in the supercritical case.

Proposition 16

The critical case is given by p

Let us plot p a and p c in Figure 4 to obtain in a different way the distinction between the supercritical case and the subcritical one. Remind we aim at estimate ψ 0,M = ψ 0 • . . . • ψ M and so we have to iterate

and ψ (2) (u) = (k + 1) log(P 0 e u + 1 -P 0 ) Heuristic : Whatever the case studied (m > 1, m < 1), we simply use two bounding functions for ψ. More precisely, for example, for u ≥ 0, until some a (intersection between the lowest homography (m < 1) and the first bissector) we bound ψ by the homographic function and then by the asymptotic direction in both cases.

Following Section 3.2.1 in a simplified version, we deduce

About random walk on the affine group and consequences

We would like to estimate

Random walk on the affine group : Consider the affine transformations

Hence the link between the random walk on the group of the homographic transformations and the one on the affine group.

Therefore, we are interested in the asymptotic properties of g 0 • g 1 • . . . • g n 's distribution which is obviously characterized by the law of ξ. As for the BPRE, we distinguish three different regimes determined by the position of E(log ξ) with respect to the origin.

For general details on random walk on affine groups the reader is referred for instance to Vervaat [START_REF] Vervaat | On a stochastic difference equation and a representation of nonnegative infinitely divisible random variables[END_REF], Brandt [START_REF] Brandt | The stochastic equation Y n+1 = A n Y n + B n with stationary coefficients[END_REF] for the one-dimensional case and to Kesten [START_REF] Kesten | Random difference equations and renewal theory for products of random matrices[END_REF], Bougerol and Picard [START_REF] Bougerol | Strict stationarity of generalized autoregressive processes[END_REF], Babillot et al. [START_REF] Babillot | The random difference equation X n = A n X n-1 + B n in the critical case[END_REF] for the d-dimensional case.

Corollary 21 Let Z L = N(0, 1) and σ 2 = var(ξ).

Confidence interval

Finally, as in the deterministic case, we deduce Theorem 5.1 For α small enough,

where

. By Jensen's inequalities,

where

Using the same kind of argument as in Corollary 6, we obtain the second part of Proposition 3.

Numerical illustration

We plot in Figure 5 (better than the first ones) of (P i ) i=1...M +1 are derived.

(3) Repeat that procedure until the budget is entirely consumed.

The goal of this algorithm is to be as close as possible to the optimal algorithm. The precise study of the proportion of the budget to use in each step shall be the purpose of a forthcoming paper and is not derive straightforward.