Charles-Edouard Bréhier
email: brehierc@cermics.enpc.fr

Tony Lelièvre
email: lelievre@cermics.enpc.fr

Mathias Rousset
email: mathias.rousset@inria.fr

ANALYSIS OF ADAPTIVE MULTILEVEL SPLITTING ALGORITHMS IN AN IDEALIZED CASE

Keywords: 1991 Mathematics Subject Classification. 65C05, 65C35, 62G30 Monte-Carlo simulation, rare events, multilevel splitting

 is a very powerful and versatile method to estimate rare events probabilities. It is an iterative procedure on an interacting particle system, where at each step, the k less well-adapted particles among n are killed while k new better adapted particles are resampled according to a conditional law. We analyze the algorithm in the idealized setting of an exact resampling and prove that the estimator of the rare event probability is unbiased whatever k. We also obtain a precise asymptotic expansion for the variance of the estimator and the cost of the algorithm in the large n limit, for a fixed k.

where a > 0 is a given threshold, such that p > 0. When a goes to infinity, the above probability goes to 0, meaning that {X ≥ a} becomes a rare event. Such problems appear in many contexts, such as molecular dynamics simulations [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF] or reliability problems with many industrial applications, for example.

Estimating a rare event using a direct Monte Carlo estimation is inefficient, as can be seen by the analysis of the relative error. Indeed, let (X i ) i∈N be a sequence of independent and identically distributed random variables with the same law as X. Then for any positive integer M , (2)

pM = 1 M M n=1
1 Xn≥a is an unbiased estimator of p: E[p M ] = p. It is also well-known that its variance is given by Var(p M ) = p(1-p) M and therefore the relative error writes:

(3)

Var(p M ) p = 1 -p M p .
Assume that the simulation of one random variable X n requires a computational cost c 0 . For a relative error of size ǫ, the cost of a direct Monte Carlo method is thus of the order (4) c 0 1p ǫ 2 p .

which is prohibitive for small probabilities (say 10 -9 ). Many algorithms devoted to the estimation of the probability of rare events have been proposed. Here we focus on the so-called Adaptive Multilevel Splitting (AMS) algorithm [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF].

Let us explain one way to understand this algorithm (see [START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF] for a more general presentation). The splitting strategy relies on the following remark. Let us introduce J intermediate levels: a 0 = 0 < a 1 < . . . < a J = a. The small probability p satisfies: p = J j=1 p j where p j = P(X > a j |X > a j-1 ), with 1 ≤ j ≤ J -1 and p J = P(X ≥ a|X > a J-1 ). In order to use this identity to build an estimator of p, one needs to (i) define appropriately the intermediate levels and (ii) find a way to sample according to the conditional distributions L(X|X > a j-1 ) to approximate each p j using independent Monte Carlo procedures.

In this article, we will be interesting in the idealized case where we assume we have a way to draw independent samples according to the conditional distributions L(X|X > a j ), j ∈ {1, . . . , J -1}. We will discuss at length this assumption below. It is then easy to check that for a given J, the variance is minimized when p 1 = . . . = p J = p 1/J , and that the associated variance is a decreasing function of J. It is thus natural to try to devise a method to find the levels a j such that p 1 = . . . = p J . In AMS algorithms, the levels are defined in an adaptive and random way in order to satisfy (up to statistical fluctuations) the equality of the factors p j . This is based on an interacting particle system approximating the quantiles P(X ≥ a) using an empirical distribution. The version of the algorithm we study depends on two parameters: n and k. The first one denotes the total number of particles. The second one denotes the number of resampled particles at each iteration: they are those with the k lowest current levels (which mean that a sorting procedure is required). Thus, the levels are defined in such a way that p j = 1 -k n and the estimator of the probability p writes:

pn,k = C n,k 1 - k n J n,k
where J n,k is the number of iterations required to reach the target level a, and C n,k ∈ 1 -k-1 n , 1 is a correction factor precisely defined below (see (10)). Notice that C n,k = 1 if k = 1.

In all the following, we will make the following assumption:

Assumption 1.1. X is a real-valued positive random variable which admits a continuous cumulative distribution function t → P(X ≤ t).

This ensures for example that (almost surely), there is no more than one particle at the same level, and thus that the resampling step in the algorithm is well defined. We will show in Section 2.3 below that this assumption can be relaxed: the continuity of the cumulative distribution function is actually only required on [0, a), in which case the AMS algorithm still yields an estimate of P(X ≥ a) (with a large inequality). From Section 3, we will always work under Assumption 1.1, and thus we will always use for simplicity strict rather than large inequalities on X (notice that under Assumption 1.1, p = P(X ≥ a) = P(X > a)).

Our aim in this article is twofold. First, we show that for any values of n and k, the estimator pn,k is an unbiased estimator of p (see Theorem 4.1): [START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF] E(p n,k ) = p.

Second, we are able to obtain an explicit asymptotic expression for the variance of the estimator pn,k in the limit of large n, for fixed k and p, and thus for the relative error (see Proposition 5.2, Equation (36)): ( 6)

Var(p n,k ) p = -log p √ n 1 + (1 -log(p)) (k -1) 2n + o 1 n 1/2
.

Thus, if we consider the cost associated to the Monte Carlo estimator based on M independent realizations of the algorithm, and if M is chosen in such a way that the relative error is of order ǫ, one ends up with the following asymptotic cost for the AMS algorithm (see Theorem 5.4, Equations ( 32) and ( 38)): for fixed k and p, in the limit of large n

(7) c 0 + c 1 log(n) ǫ 2 (log(p)) 2 -log(p) + 1 n -log(p) (k -1) + 1 2 (log(p)) 2 - 1 2 (log(p)) 3 + o 1 n
where c 0 denotes the cost for drawing one sample according to the conditional distributions L(X|X > x) (assumed to be independent of x, to simplify), and c 1 log(n) is the cost associated with the sorting procedures involved in the algorithm. Here again the Landau symbol o in the above depend non-uniformly on k and p. The two results ( 6) and [START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] should be compared to the corresponding formulae for direct Monte Carlo simulation (3) and (4) above. From these results, we conclude that, in this asymptotic regime: (i) the choice k = 1 is the optimal one in terms of variance and cost (ii) AMS yields better result than direct Monte Carlo if

1 + c 1 c 0 log(n) (log p) 2 -log p < 1 -p p
which will be the case for sufficiently small probability p. The assumption that we are able to sample according to the conditional distributions L(X|X > x) (idealized setting) is a severe limitation, from a practical viewpoint. Let us make three comments about this idealized setting. First, to the best of our knowledge, all the theoretical results which have been obtained so far in the literature [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF][START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF][START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] rely on such an assumption. It is believed that the qualitative conclusions obtained under this assumption are still meaningful for the actual algorithm used in practice, where the conditional distributions L(X|X > a j-1 ) are only approximately sampled (using for example Metropolis Hastings procedures). In some sense, in the idealized setting, one studies the optimal performance one could reach with this algorithm. Second, we will describe in Section 2.2 situations where this assumption makes sense in practice: this is in particular the case in the 1d dynamic setting described in Section 2.2.2. Third, in a paper in preparation [START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF], we will actually show that it is possible to obtain an unbiased estimator of p using AMS in a very general setting. In other words, the idealized setting is crucial to obtain estimates on the variances and the costs of the algorithm, but not to prove the unbiasedness property for the estimator of the rare event probability.

Let us now review the results known from the literature on AMS. As mentioned above, the AMS algorithm has been introduced in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF], where it is proven that their estimator pn,k indeed converges almost surely to p (in the limit n → ∞). A central limit theorem is also provided, in the asymptotic regime where k n = p 0 is fixed. In [START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF], the authors consider the case k = 1, prove the unbiasedness of the estimator and analyze the variance. In [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF], the authors analyze the large n limit, with fixed k n = p 0 . In summary, our results differ from what has been proven before in two ways: we prove that pn,k is an unbiased estimator for any k ≥ 1 (and not only k = 1) and we analyze the variance and the cost at a fixed k (in the limit n → ∞) and show that in this regime, k = 1 is optimal.

In addition to the new results presented in this paper, we would like to stress that the techniques of proof we use seem to be original. The main idea is to consider (under Assumption 1.1) the family (P (x)) x∈[0,a] of conditional probabilities:

(8) P (x) = P(X > a|X > x),
and to define associated estimators pn,k (x) thanks to the AMS algorithm. We are then able to derive an explicit functional equation on E[p n,k (x)] as a function of x; which follows from the existence of an explicit expression for the distribution of order statistics of independent random variables. We can then check the unbiased property (see Theorem 4.1) E[p n,k (x)] = P (x). Finally noting that P (0) = p, we get [START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF]. To analyze the computational cost (Theorem 5.4), we follow a similar strategy, with several more technical steps. First, we derive functional equations for the variance and the mean number of iterations in the algorithm. In general we are not able to give explicit expressions to the solutions of these equations, and we require the following auxiliary arguments:

(i) We show how one can relate the general case to the so-called exponential case, when X has an exponential distribution with mean 1; (ii) We then prove that the solutions of the functional equations are solutions of linear Ordinary Differential Equations of order k; (iii) We finally get asymptotic results on the solutions to these ODEs in the limit n → +∞, with fixed values of p and k. The paper is organized as follows. In Section 2, we introduce the AMS algorithm and discuss the sampling of the conditional distributions L(X|X > z). In Section 3, we show how to relate the case of a general distribution for X to the case when X is exponentially distributed. We then prove that the algorithm is well defined, in the sense that it terminates in a finite number of iterations (almost surely). In Section 4, we show one of our main result, namely Theorem 4.1 which states that the estimators pn,k (x) are unbiased. Finally, in Section 5, we study the cost of the algorithm, with asymptotic expansions in the regime where p and k are fixed and n goes to infinity. The proofs of the results of Section 5 are postponed to Section 6.

The Adaptive Multilevel Splitting algorithm

After presenting the AMS algorithm in Section 2.1, we discuss the fundamental assumption that we know how to sample according to the conditional distributions L(X|X > z) in Section 2.2. We will in particular show that this assumption is actually practical at least in one setting: the one-dimensional dynamic setting presented in Section 2.2.2. Finally, we discuss Assumption 1.1 in Section 2.3 and show that it can be replaced by a less stringent hypothesis. This is particularly useful in the framework of the high-dimensional dynamic setting described in Section 2.2.3.

We would like to stress that Sections 2.2 and 2.3 discuss practical aspects of the implementation of AMS and can be skipped if the reader is only interested in the two main results about the unbiasedness and cost of AMS, presented in Sections 4 and 5.

2.1. Description of the algorithm. We fix a total number n of particles, as well as 1 ∈ {1, . . . , n -1} the number of resampled particles at each iteration of the algorithm. In the sequel, when we consider a random variable X j i , the subscript i denotes the index in {1, . . . , n} of a particle, and the superscript j denotes the iteration of the algorithm.

In the algorithm below and in the following, we use classical notations for k-th order statistics. For Y = (Y 1 , . . . , Y n ) an ensemble of independent and identically distributed (i.i.d.) real valued random variables with continuous cumulative distributions function, there exists almost surely a unique (random) permutation σ of {1, . . . , n} such that Y σ(1) < . . . < Y σ(n) . For any k ∈ {1, . . . , n}, we then use the classical notation Y (k) = Y σ(k) to denote the k-th order statistics of the sample Y .

For any x ∈ [0, a], we define the Adaptive Multilevel Splitting algorithm as follows (in order to approximate p = P(X ≥ a) one should take x = 0, but we consider the general case x ∈ [0, a] for theoretical purposes).

Algorithm 2.1 (Adaptive Multilevel Splitting).

Initialization: Define Z 0 = x. Sample n i.i.d. realizations X 0 1 , . . . , X 0 n , with the law L(X|X > x). Define Z1 = X 0 (k) , the k-th order statistics of the sample X 0 = (X 0 1 , . . . , X 0 n ), and σ 1 the (a.s.) unique associated permutation:

X 0 σ 1 (1) < . . . < X 0 σ 1 (n) . Set j = 1.
Iterations (on j ≥ 1): While Z j < a:

Conditionally on Z j , sample k new independent 1 random variables (χ j 1 , . . . , χ j k ), according to the law L(X|X > Z j ). Set

X j i = χ j (σ j ) -1 (i) if (σ j ) -1 (i) ≤ k X j-1 i if (σ j ) -1 (i) > k.
In other words, the particle with index i is killed and resampled according to the law L(X|X > Z j ) if X j-1 i ≤ Z j , and remains unchanged if X j-1 i > Z j . Notice that the condition (σ j ) -1 (i) ≤ k is equivalent to i ∈ σ j (1), . . . , σ j (k) .

Define Z j+1 = X j (k) , the k-th order statistics of the sample X j = (X j 1 , . . . , X j n ), and σ j+1 the (a.s.) unique2 associated permutation: X j σ j+1 (1) < . . . < X j σ j+1 (n) . Finally increment j ← j + 1. End of the algorithm: Define J n,k (x) = j -1 as the (random) number of iterations. Notice that J n,k (x) is such that Z J n,k (x) < a and Z J n,k (x)+1 ≥ a.

The estimator of the probability p x is defined by

(9) pn,k (x) = C n,k (x) 1 - k n J n,k (x)
, with

C n,k (x) = 1 n Card i; X J n,k (x) i ≥ a . (10)
Notice that C n,1 (x) = 1. More generally, C n,k (x) ∈ n-k+1 n , . . . , n-k+k n . Since we are interested in the algorithm starting at x = 0, we introduce the notation (11) pn,k = pn,k (0).

We finally stress that the computations of the sampled random variables (X 0 i ) 1≤i≤n for the initialization and of the (χ j i ) 1≤i≤k for each iteration j can be made in parallel. 2.2. Sampling from the conditional distributions L(X|X > z). At each iteration of the algorithm, we need to sample k random variables according to conditional distributions L(X|X > z), with z taking values in the sequence (Z j ) 0≤j≤J n,k (x) . As explained above, we develop our theoretical analysis of the properties of the algorithm (bias, variance and computational cost) in the idealized situation where it is possible to sample according to these conditional distributions L(X|X > z) for any z ∈ [0, a]. From a practical point of view, this assumption is generally unrealistic. One possible situation where it is realistic is the dynamic setting presented in the present section (in contrast with the static setting).

2.2.1.

The static setting and the exponential case. In a general framework, there is no simple way to sample the distributions L(X|X > z). In practice, this can be done thanks to a Metropolis-Hastings procedure, see for example [START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF]. Of course, this introduces a bias and correlations between the particles at each iteration (compared with the idealized algorithm studied in this paper). This bias and these correlations asymptotically vanish when the number of iterations in the Metropolis-Hastings procedure goes to infinity. The error analysis associated to this procedure is out of the scope of this paper. There is a simple example where it is actually possible to sample the distributions L(X|X > z), namely if X is exponentially distributed. Indeed, if X has exponential law E(1) with mean 1, then the conditional distribution L(X|X > x) = L(X + x) is a shifted exponential variable, for any x > 0. In the following, we will refer to this situation as the exponential case. Of course, this has no practical interest since in this case, p = P(X ≥ a) = exp(-a) is analytically known. However, this particular case plays a crucial role in the analysis hereafter, since as will be precisely explained in Section 3.3, the study of the general case can be reduced to the study of the exponential case after some change of variable. This trick was already used in the original papers [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF][START_REF] Cérou | Adaptive particle techniques and rare event estimation[END_REF].

2.2.2. Dynamic setting in dimension 1. In the one-dimensional dynamic setting, X is defined as

X = sup 0≤t≤τ Y t
where (Y t ) 0≤t≤τ is a strongly Markovian time-homogeneous random process with values in R, and τ is a stopping time. In this setting, the conditional distribution L(X|X > x) is easily sampled: it is the law of sup 0≤t≤τ Y x t , where Y x t denotes the stochastic process (Y t ) t≥0 which is such that Y 0 = x. Having in mind applications in molecular dynamics [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF], a typical example is when (Y t ) t≥0 satisfies a Stochastic Differential Equation

dY x t = f (Y x t )dt + 2β -1 dB t , Y x 0 =
x, with smooth drift coefficient f and inverse temperature β > 0. The stopping time is for example [START_REF] Asmussen | Stochastic simulation: Algorithms and analysis Springer[END_REF], and for some ǫ > 0 and one can then consider X = sup 0≤t≤τ 0 Y 0 t . Let us consider the target level a = 1. The AMS algorithm then yields an estimate of P(X ≥ 1) the probability that the stochastic process starting from 0 reaches the level 1 before the level -ǫ. Such computations are crucial to compute transition rates and study the so-called reactive paths in the context of molecular dynamics, see [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF].

τ x = inf {t ≥ 0; Y x t < -ǫ or Y x t > 1 + ǫ} , for x ∈ [0,
Notice that in practice, a discretization scheme must be employed, which makes the sampling of the conditional probabilities more complicated. Another point of view on the AMS algorithm is then required in order to prove the unbiasedness of the estimator of the probability p, see [START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF].

Remark 2.2. The exponential case can be obtained from a dynamic setting. Indeed, consider the following stochastic process: a particle starts at a given position x, moves on the real line with speed p = +1 on the random interval [0, τ ] where τ ∼ E(1) is exponentially distributed, and with speed p = -1 on the interval (τ, +∞). More precisely,

p x t = +1 for 0 ≤ t ≤ τ -1 for t > τ q x t = x + t for 0 ≤ t ≤ τ x + τ -t for t > τ.
Notice that (p x t , q x t ) t≥0 is a Markov process such that (q x t ) t≥0 is continuous. Then for any initial condition x and any given threshold a > x we have P(sup t≥0 q x t > a) = P(τ > ax) = exp(xa). In particular, X = sup t≥0 q 0 t = sup 0≤t≤τ q 0 t is an exponential random variable with parameter 1. 2.2.3. Dynamic setting in higher dimension. Let us consider again a strongly Markovian time-homogeneous stochastic process (Y t ) t≥0 , but with values in R d for d ≥ 2. In this case, the levels need to be defined using a (continuous) function ξ : R d → R, sometimes called a reaction coordinate in the context of molecular dynamics.

Let us focus for simplicity on the case when (Y t ) t≥0 is solution of the stochastic differential equation (SDE):

(12)

dY x t = -∇V (Y x t )dt + 2β -1 dW t , Y x 0 =
x, with smooth potential V , inverse temperature β > 0 and (W t ) t≥0 a d-dimensional Wiener process. Let us consider two disjoint closed subsets A and B of R d . Let us define the stopping time

τ x = min(τ x A , τ x B ) where (13) τ x A = inf {t ≥ 0; Y x t ∈ A} and τ x B = inf {t ≥ 0; Y x t ∈ B} .
Let us assume that the function ξ is such that

A = {x; ξ(x) ≤ 0} and B = {x; ξ(x) ≥ 1}.
We then set, for a fixed initial condition

x 0 ∈ R d \ (A ∪ B), (14)
X = sup 0≤t≤τ ξ(Y x0 t ).
Let us set a = 1 as the target level. In this case, the probability p = P(X ≥ 1) = P(τ x0 B < τ x0 A ) is the probability that the path starting from x 0 reaches B before A. As explained above, this is a problem of interest in molecular dynamics for example, to study reactive paths and compute transition rates in high dimension, typically when A and B are metastable regions for (Y t ) t≥0 .

The problem to apply the AMS algorithm is again to sample according to the conditional distributions L(X|X > z). A natural idea is to use the following branching procedure in the resampling step at the j-th iteration: to build one of the new k trajectories, one of the (nk) remaining trajectories is chosen at random, copied up to the first time it reaches the level {x; ξ(x) = Z j } and then completed independently from the past up to the stopping time τ . The problem is that this yields in general a new trajectory which is correlated to the copied one through the initial condition on the level set {x; ξ(x) = Z j }. Indeed, in general, given

x 1 = x 2 such that ξ(x 1 ) = ξ(x 2 ), the laws of sup 0≤t≤τ ξ(Y x1 t ) and sup 0≤t≤τ ξ(Y x2 t )
are not the same. As a consequence, it is unclear how to sample L(X|X > z), except if we would be able to build a function ξ such that the law of sup 0≤t≤τ ξ(Y x t ) only depends on ξ(x). This is actually the case if ξ is the so-called committor function associated to the dynamics (12) and the two sets A and B.

(15)      -∇V • ∇ξ + β -1 ∆ξ = 0 in R d \ (A ∪ B), ξ(x) = 0 for x ∈ A, ξ(x) = 1 for x ∈ B.
Proposition 2.4. Assume that ξ is the committor function, solution of (15). For any

x ∈ R d \ (A ∪ B), we set X x = sup 0≤t≤τ ξ(Y x t ). Notice that X x is a random variable with values in [ξ(x), 1]. We have P(X x > z) = ξ(x) z for any x ∈ R d \ (A ∪ B) and z ∈ [ξ(x), 1). In particular, if ξ(x 1 ) = ξ(x 2 ), then we have the equality L(X x1 ) = L(X x2 ). Moreover, for any x 1 , x 2 with 0 < ξ(x 1 ) ≤ ξ(x 2 ) < 1, L(X x1 |X x1 > ξ(x 2 )) = L(X x2 ).
This previous Proposition (which is proven below) fully justifies the branching procedure described above to sample L(X x0 |X x0 > Z j ) at the j-th iteration of the algorithm: pick at random one of the (nk) remaining trajectories (say (Y t ) t≥0 ), copy it up to the first time it reaches the level {x; ξ(x) = Z j } (let us denote Y σ the first hitting point of this level) and then complete the trajectory, independently from the past. By the strong Markov property, this yields a new X sampled according to L(X Yσ ) which is indeed

L(X x0 |X x0 > Z j ), since ξ(Y σ ) = Z j .
Remark 2.5. As already mentioned in the previous Section, in practice, the SDE (12) is discretized in time, say with a timestep ∆t. Then, the branching procedure consists in copying the selected trajectory up to the first time index n such that ξ(Y n∆t ) > z, and then to complete it independently from the past. This introduces a difference compared to the continuous in time situation considered above, since doing so we do not sample according to the conditional distribution L(X|X > z). To treat this situation, one needs to resort to other techniques to analyze the algorithm, see [START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF]. In particular, one can show that the algorithm still yields an unbiased estimator of p, using very different techniques of proof than those presented in this paper. This approach is also useful to treat non-homogeneous in time Markov processes.

At this stage, we can thus conclude that in the high-dimensional dynamic setting, if the committor function is known, the AMS algorithm can be practically implemented, and that it enters the framework of this paper. There are however two difficulties, that we will now discuss.

First, the random variable X = sup 0≤t≤τ ξ(Y x0 t ), where ξ is the committor function, does not satisfy Assumption 1.1: we have p = P(X ≥ 1) = P(X = 1) > 0 (X takes values in [0, 1]) and therefore, the cumulative distribution t → P(X ≤ t) is not continuous at t = 1. More precisely, from Proposition 2.4, we have

: ∀t ∈ [0, 1), P(X ≤ t) = 1 -ξ(x0) t +
and ∀t ≥ 1, P(X ≤ t) = 1. This is actually not a problem, as explained in the next Section in a general setting: the continuity of the cumulative distribution function is only required over [0, 1) (or more generally over [0, a) in the general case when the target level is a).

The second difficulty is that knowing the committor function is actually a very strong assumption. Computing ξ solution to (15) is actually impossible in practice since this is a high-dimensional PDE. Moreover, if ξ was known, then we would actually know the small probability we want to estimate since p = P(X ≥ 1) = ξ(x 0 ). This is a consequence of the well-known probabilistic representation to solutions to (15): Proposition 2.6. Recall the definitions (13) of the stopping times τ x A and τ x B . Then, if ξ is the committor function associated to the SDE (12) and the sets A and B, then, for any

x ∈ R d (16) ξ(x) = P(τ x B < τ x A
). Thus, this high-dimensional dynamic case with known committor function should also be considered as an idealized setting, which is only useful for theoretical purposes, in order to study the best performance we could expect for the AMS algorithm.

We end up this Section with a proof of Proposition 2.4. Proof of Proposition 2.4: Let us consider ξ satisfying (15), Y x t solution to (12) and X x = sup 0≤t≤τ ξ(Y x t ). For any given z ∈ (0, 1), and any x ∈ R d , let us introduce

τ x z = inf {t ≥ 0; ξ(Y x t ) ≥ z} . One easily checks the identity P(X x ≥ z) = P(X x > z) = P(τ x z < τ x A )
. By continuity of ξ and of the trajectories of the stochastic process (Y t ) t≥0 , and by the strong Markov property at the stopping time τ x z , we get for any x and any z ∈ (ξ(x), 1)

ξ(x) = P(τ x B < τ x A ) = E 1 τ x z <τ x A 1 τ x B <τ x A = E 1 τ x z <τ x A 1 τ Y x τ x z B <τ Y x τ x z A = E 1 τ x z <τ x A ξ(Y x τ x z ) = P(τ x z < τ x A )z.
This identity proves the first claim of the Proposition. Moreover, since the law of Y x depends on x only through ξ(x), we also more generally get

L(X x1 |X x1 > ξ(x 2 )) = L(X x2 ) for any x 1 , x 2 with 0 < ξ(x 1 ) ≤ ξ(x 2 ) < 1,.
2.3. About Assumption 1.1. In this section, we show that Assumption 1.1 is actually too stringent. If one assumes the following Assumption 2.7. X is a real-valued positive random variable such that t ∈ [0, a) → P(X ≤ t) is continuous, then the Algorithm 2.1 is well defined, and all the results presented below hold. In particular, the estimator pn,k is an unbiased estimator of p = P(X ≥ a). We notice that Assumption 2.7 is indeed more natural than Assumption 1.1 since the AMS algorithm only applies a resampling procedure with conditional distributions L(X|X > Z j ) to realizations such that Z j ∈ [0, a): this is why the continuity of the cumulative distribution function t → P(X ≤ t) is actually only required over [0, a).

The argument to show that one can recover the setting of Assumption 1.1 assuming only Assumption 2.7 is the following coupling argument. In Lemma 2.8 below, it is proven that there exists a random variable X such that:

• X satisfies Assumption 1.1;

• for any z ∈ [0, a), X ≤ z is equivalent to X ≤ z and L( X| X > z) = L(X|X > z);
• X ≥ a is equivalent to X ≥ a and thus, in particular, P( X ≥ a) = P(X ≥ a) = p, the probability to be estimated. The last two properties show that running the AMS algorithm on X is equivalent to running the AMS algorithm on X: the iterations, the stopping criterion and the estimator are the same. The theory developed in this paper (unbiased estimator, analysis of the cost and of the computational cost) is then applied to the algorithm applied to the auxiliary random variable X instead of X, which is completely equivalent to the algorithm applied to X.

In all the following, for simplicity, we will always assume that Assumption 1.1 holds, keeping in mind that it can be relaxed to Assumption 2.7. Thus, inequalities which involve the random variable X can be changed from large to strict without modifying the associated events (almost surely).

We end this Section with a Lemma which defines the random variable X as a function of X.

Lemma 2.8. Let X be a random variable satisfying Assumption 2.7, and let us define

X = X1 X<a + a U 1 X≥a ,
where U is a random variable independent of X and uniformly distributed on (0, 1). Then, (i) X satisfies Assumption 1.1, (ii) for any z ∈ [0, a), X ≤ z is equivalent to X ≤ z, and the two laws L( X| X > z) and L(X|X > z) coincide on (z, a) and (iii) X ≥ a is equivalent to X ≥ a.

Proof of Lemma 2.8: Since a/U > a, it is easy to check the items (ii) and (iii). Let us now consider the cumulative distribution of X.

For t < a, P( X ≤ t) = P(X ≤ t) and thus t → P( X ≤ t) is continuous for t ∈ [0, a) by Assumption 2.7. For t > a, P(

X ≤ t) = P(X < a) + P(U ≥ a/t, X ≥ a) = P(X < a) + P(X ≥ a) 1 -a t = 1 -a t P(X ≥ a) and thus t → P( X ≤ t) is continuous for t ∈ (a, +∞).
Finally, with these expressions one easily checks left and right continuity of t → P( X ≤ t) at a. This concludes the proof of the fact that X satisfies Assumption 1.1, and thus the proof of Lemma 2.8.

Reduction to the exponential case and well-posedness of the algorithm

The aim of this Section (see Section 3.3) is to prove the well-posedness of the algorithm, namely the fact that J n,k (x) is almost surely finite, when the probability p = P (0) to estimate is positive. The argument is based on the fact that the general case is related to the exponential case through a change of variable, see Section 3.2 (this will be instrumental in the rest of the paper). Section 3.1 first gives a few notation that will be useful below.

Notation.

3.1.1. General notation. We will use the following set of notations, associated to the random variable X satisfying Assumption 1.1.

We denote by F the cumulative distribution function of the random variable X: F (t) = P(X ≤ t) for any t ∈ R, and F (0) = 0. From Assumption 1.1, the function F is continuous. Notice that it ensures that if Y is an independent copy of X, then P(X = Y ) = 0, and thus, in the algorithm, there is only at most one sample at a given level.

We recall that our aim is to estimate the probability

(17) p = P(X ≥ a) = P(X > a),
given a threshold a > 0. More generally, we define for any x ∈ [0, a] (18) P (x) = P(X > a|X > x), so that we have p = P (0). Notice that P (a) = 1.

For any x ∈ [0, a], L(X|X > x) admits a cumulative distribution function F (.; x), which satisfies: for any y ∈ R,

(19) F (y; x) = F (y) -F (x) 1 -F (x) 1 y≥x .
While the variable x always denotes the parameter in the conditional distributions L(X|X > x), we use the variable y as a dummy variable in the associated densities and cumulative distribution functions. By Assumption 1.1, 1 y≥x in the definition above can be replaced with 1 y>x . Notice that F (y; 0) = F (y). Moreover, with these notations, we have (20)

P (x) = 1 -F (a; x).
An important tool in the following is the family of functions: for any x ∈ [0, a] and any y ∈ R

(21) Λ(y; x) = -log(1 -F (y; x)) ∈ [0, +∞] Λ(y) = Λ(y; 0) = -log(1 -F (y)).
We remark the following identity: for 0 ≤ x ≤ y ≤ a,

Λ(y; x) = log(1 -F (x)) -log(1 -F (y)) = Λ(y) -Λ(x).
3.1.2. Specific notation when X admits a density f . In some places, we will assume that X admits a density f with respect to the Lebesgue measure (which indeed implies Assumption 1.1). This assumption is in particular satisfied in the exponential case (namely when X is exponentially distributed), which we will consider in several arguments below to study the bias and the computational cost. If X admits a density f , the law L(X|X > x) of X conditionally on {X > x} also admits a density f (.; x), which satisfies: for any y ≥ 0,

(22) f (y; x) = f (y) 1 -F (x) 1 y≥x .
Notice that f (y) = f (y; 0). We finally introduce some notations about order statistics of samples of i.i.d. real random variables. If X admits a density f , then the k-th order statistics of an i.i.d. sample (X 1 , . . . , X n ) (distributed according to the law of X) admits a density f n,k which satisfies: for any y ∈ R,

f n,k (y) = k n k F (y) k-1 f (y) 1 -F (y) n-k .
The associated cumulative distribution function is

F n,k (y) = y -∞ f n,k (z)dz.
Likewise, we introduce notations for the density and the cumulative distribution function of order statistics for the law L(X|X > x): when 0 ≤ x ≤ y ≤ a we set (23)

f n,k (y; x) = k n k F (y; x) k-1 f (y; x) 1 -F (y; x) n-k , F n,k (y; x) = y x f n,k (z; x)dz.
3.2. Reduction to the exponential case. One of the key tool in the following is the reduction of the general case to the exponential case, thanks to the use of the function Λ, defined by (21). We recall that the exponential case refers to the case when X is distributed according to the exponential law with parameter 1, see Section 2. Let us first state a result on the algorithm without any stopping criterion.

Proposition 3.2. Let us consider the sequence of random variables ((X j i ) 1≤i≤n , Z j ) j≥0 generated by the AMS algorithm 2.1 without any stopping criterium.Set Y j i = Λ(X j i ) and S j = Λ(Z j ). Then we have the following properties.

(i) For any j ≥ 0, (Y j i -S j ) 1≤i≤n is a family of i.i.d. exponentially distributed random variables, with parameter 1.

(ii) For any j ≥ 1, (Y j i -S j ) 1≤i≤n is independent of (S l -S l-1 ) 1≤l≤j . (iii) The sequence (S j -S j-1 ) j≥1 is i.i.d.. As a consequence, in law, the sequence ((Y j i ) 1≤i≤n , S j ) j≥0 is equal to the sequence of random variables ((X j i ) 1≤i≤n , Z j ) j≥0 obtained by the realization of the AMS algorithm without any stopping criterion in the exponential case, with initial condition Z 0 = Λ(x).

Proof of Proposition 3.2: The last assertion is a direct consequence of the three former items and of the fact that in the exponential case the function Λ is the identity mapping.

Item (iii) is a consequence of items (i) and (ii). It remains to prove jointly those two items, which is done by induction on j. When j = 0, the result follows from the way the algorithm is initialized: for each

1 ≤ i ≤ n, Y 0 i = Λ(X 0 i
) is exponentially distributed (see Lemma 3.1), and the independence property is clear. Assuming that the properties (i) and (ii) are satisfied for all k ≤ j, it is sufficient to prove that:

(i) (Y j+1 i -S j+1 ) 1≤i≤n are i.i.d. and exponentially distributed with mean 1;

(ii) (Y j+1 i -S j+1 ) 1≤i≤n is independent of (S l -S l-1 ) 1≤l≤j+1 .
For any positive real numbers y 1 , . . . , y n , and any s 1 , . . . , s j+1 , this is equivalent to proving that (24)

A : = P Y j+1 1 -S j+1 > y 1 , . . . , Y j+1 n -S j+1 > y n , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 = exp (-(y 1 + • • • + y n )) P(S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 ).
We want to decompose the probability with respect to the value of σ j+1 (i)) 1≤i≤k . We recall that almost surely we have

Y j σ j+1 (1) < . . . < Y j σ j+1 (k) = S j+1 < Y j σ j+1 (k+1) < . . . < Y j σ j+1 (n) .
In fact, in order to preserve symmetry inside the groups of resampled and not-resampled particles, we decompose over the possible values for the random set σ j+1 (1) . . . σ j+1 (k) . We thus compute a sum over all partitions {1, . .

P Y j+1 i -S j+1 > y i ; i ∈ I + , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 i∈I- P Λ(χ j+1 ) -Λ(Z j+1 ) > y i .
In the last line, we used the fact that by construction of the algorithm, on the event we consider, namely σ j+1 (1), . . . , σ j+1 (k) = I -, we have the equality of random variables

Y j+1 i = Λ χ j+1 (σ j+1 ) -1 (i) , i ∈ I -.
The latter are i.i.d. and independent of all the other random variables used at this stage of the algorithm. Moreover:

i∈I- P Λ(χ j+1 ) -Λ(Z j+1 ) > y i = i∈I- exp -y i .
Let us now introduce a notation: for

I -⊂ {1, • • • n} we set M j I-:= max Y j i -S j ; i ∈ I -.
Remark that on the event σ j+1 (1), . . . , σ j+1 (k) = I -, almost surely, we have Y j i = Y j+1 i and S j+1 -S j = M j I-. Using the independence properties of the Y j i -S j , 1 ≤ i ≤ n (and thus of M j I-) from the induction hypothesis (ii), we obtain

P Y j+1 -S j+1 > y i ; i ∈ I + , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 = P Y j i -S j -M j I-> y i ; i ∈ I + , S 1 -S 0 > s 1 , . . . , M j I-> s j+1 = P Y j i -S j -M j I-> y i ; i ∈ I + , M j I-> s j+1 P S 1 -S 0 > s 1 , . . . , S j -S j-1 > s j .
We can then integrate the Y j i -S j , i ∈ I + using the induction hypothesis (i):

P Y j i -S j -M j I-> y i ; i ∈ I + , M j I-> s j+1 = E i∈I+ exp -y i -M j I-1 M j I - >s j+1 = i∈I+ exp -y i E exp -(n -k)M j I-1 M j I - >s j+1
The proof is now complete since:

A = P Y j+1 1 -S j+1 > y 1 , . . . , Y j+1 n -S j+1 > y n , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 = I -⊂{1,...,n} Card(I-)=k i∈I- exp -y i i / ∈I- exp -y i E exp -(n -k)M j I-1 M j I - >s j+1 P S 1 -S 0 > s 1 , . . . , S j -S j-1 > s j = exp n i=1 y i P S 1 -S 0 > s 1 , . . . , S j -S j-1 > s j I -⊂{1,...,n} Card(I-)=k E exp -(n -k)M j I-1 M j I - >s j+1 ,
which proves (24).

In particular, taking y 1 = . . . = y n = s 1 = . . . = s j = 0, we see that

I -⊂{1,...,n} Card(I-)=k E exp -(n -k)M j I-1 M j I - >s j+1 = P(S j+1 -S j > s j+1 ).
This concludes the proof of Proposition 3.2.

The next Lemma shows that the AMS algorithm applied to X with target level a and the AMS algorithm applied to Λ(X) with target level Λ(a) stop at the same iteration. Lemma 3.3. Let us consider the sequence of random variables ((X j i ) 1≤i≤n , Z j ) j≥0 generated by the AMS algorithm 2.1 without any stopping criterium, and set Y j i = Λ(X j i ) and S j = Λ(Z j ). For any α > 0, almost surely, S j ≥ Λ(α) = Z j ≥ α .

Proof of Lemma 3.3: : First, Λ is non-decreasing so that Z j ≥ α ⊂ S j ≥ Λ(α) . Moreover, one easily checks that S j ≥ Λ(α) ∩ Z j < α ⊂ S j = Λ(α) . From Proposition 3.2 we know that S j = S 0 + j ℓ=1 S ℓ -S ℓ-1 admits a density with respect to the Lebesgue measure, since the S ℓ -S ℓ-1 have the density of the k-th order statistics of independent exponentially distributed random variables with parameter 1. Therefore P S j > Λ(α) = Z j > α = 0.

A direct corollary of Proposition 3.2 and Lemma 3.3 is that the original problem reduces to the exponential case.

Corollary 3.4. Consider the sequences of random variables (X j i ) 1≤i≤n,0≤j≤J n,k (x) and (Z j ) 0≤j≤J n,k (x)+1 generated by the AMS algorithm 2.1. Set Y j i = Λ(X j i ) and S j = Λ(Z j ).

Then, the sequences (Y j i ) 1≤i≤n,0≤j≤J n,k (x) and (S j ) 0≤j≤J n,k (x)+1 are equal in law to the sequences of random variables (X j i ) 1≤i≤J n,k (Λ(x)) and (Z j ) 0≤j≤J n,k (Λ(x))+1 obtained by the realization of the AMS algorithm in the exponential case, with initial condition Z 0 = Λ(x) and target level Λ(a).

Finally, in the next Sections, we need the following result, which is a consequence of Proposition 3.2. Corollary 3.5. For any j ≥ 0, conditionally on Z j , the random variables (X j i ) 1≤i≤n are i.i.d. with law L(X|X > Z j ).

Proof: Thanks to Proposition 3.2, we see that (Λ(X j i )-Λ(Z j )) 1≤i≤n are i.i.d. and exponentially distributed with mean 1. Since Λ(X j i ) -Λ(Z j ) = Λ(X j i ; Z j ), we observe that for any x 1 , . . . , x n ∈ [Z j , +∞),

P(X j 1 > x 1 , . . . , X j n > x n |Z j ) = P(Λ(X j 1 ) > Λ(x 1 ), . . . , Λ(X j n ) > Λ(x n )|Λ(Z j )) = P(Λ(X j 1 ) -Λ(Z j ) > Λ(x 1 ) -Λ(Z j ), . . . , Λ(X j n ) -Λ(Z j ) > Λ(x n ) -Λ(Z j )|Λ(Z j )) = exp(-(Λ(x 1 ) -Λ(Z j )) . . . exp(-(Λ(x n ) -Λ(Z j )) = (1 -F (x 1 ; Z j )) . . . (1 -F (x n ; Z j )).
This concludes the proof.

3.3.

Well-posedness of the algorithm. To ensure that the algorithm giving an estimator of the probability is well-defined, namely that it gives a result after a finite number of steps, we prove in this Section that J n,k (x) is almost surely finite, when the probability p = P (0) is positive. The proof relies on the reduction to the exponential case explained in the previous Section. Proposition 3.6. Suppose p = P (0) > 0. Then for any x ∈ [0, a], we have P (x) = P(X > a|X > x) > 0, and for any integers n and k with 1 ≤ k < n, the number of iterations in the AMS algorithm is almost surely finite: J n,k (x) < +∞ a.s.

Proof of Proposition 3.6:

To prove this result, we consider the AMS Algorithm 2.1 without any stopping criterion (namely the condition Z j ≥ a). As a consequence, we define sequences of random variables with the iteration index j ∈ N: we get (X j i ) j≥0 , for any i ∈ {1, . . . , n} and (Z j ) j≥0 . Proposition 3.6 is then equivalent to the following statement: almost surely, j ≥ 0; Z j ≥ a = ∅.

Thanks to Proposition 3.2, we write for any j ≥ 0

S j = S 0 + j ℓ=1 R ℓ ,
where R ℓ = S ℓ -S ℓ-1 are independent and identically distributed positive random variables, satisfying ER ℓ ∈ (0; +∞). Indeed,

ER 1 ≤ E max 1≤i≤n Y 0 i ≤ E 1≤i≤n Y 0 i = n,
where we recall that Y 0 i = Λ(X 0 i ) are independent and exponentially distributed with parameter 1. To prove that ER 1 > 0, we write

ER 1 ≥ E min 1≤i≤n Y 0 i = 1/n,
since it is easily checked that min 1≤i≤n Y 0 i has an exponential distribution, with mean 1/n. By the Strong Law of Large Numbers, when j → +∞, we have the almost sure convergence

S j j → ER 1 ,
which yields S j → +∞, almost surely, when j → +∞. As a consequence, almost surely, there exists some j ∈ N such that S j ≥ Λ(a). Using Lemma 3.3, this then implies that Z j ≥ a, and that J n,k (x) < +∞.

In the case k = 1, following the ideas in the proof of Proposition 3.6, one can easily identify the law of the number of iterations (see [START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] for a similar result).

Proposition 3.7. The random variable J n,1 (x) has a Poisson distribution with mean -n log(P (x)).

Proof of Proposition 3.7: In the case k = 1, R 1 = min 1≤i≤n Y 0 i has an exponential distribution with mean 1/n. We recall that (R ℓ = S ℓ -S ℓ-1 ) ℓ≥0 is a sequence of independent and identically distributed random variables.

Let us introduce the Poisson process, with intensity n, associated with the sequence of independent and exponentially distributed increments (R j ) ℓ≥0 :

P t = +∞ ℓ=0 1 S ℓ ≤t .
Since S 0 = Λ(x), we identify that J n,1 (x) = P Λ(a)-Λ(x) = P -log(P (x)) .

The result follows since for any t ≥ 0 P t has a Poisson distribution with mean nt.

The estimator pn,k (x) is unbiased

Recall that X satisfies Assumption 1.1. Let us fix a total number of replicas n, as well as k ∈ {1, . . . , n -1} the number of killed and resampled replicas at each iteration. Given a level a > 0, we recall that the estimator of the conditional probability P (x) = P(X > a|X > x) for each value of x ∈ [0, a] is pn,k (x), defined by (9). We introduce the following notation:

(25) p n,k (x) = E[p n,k (x)].
Recall that we are specifically interested in estimating the probability p = P (0), and that the introduction of P (x) for x in the interval [0, a] is a tool to prove that the estimator is unbiased. We write the result in its full generality, and then specify it to the estimation of p. Theorem 4.1. For any k ∈ {1, . . . , n -1}, for any a > 0, such that p = P(X > a) = P (0) > 0, and any x ∈ [0, a], pn,k (x) is an unbiased estimator of the conditional probability P (x):

(26) E[p n,k (x)] = P (x).
In particular, when x = 0, we have E[p n,k ] = p.

From Section 3.2 (see Corollary 3.4), it is sufficient to prove the result in the exponential case. Indeed, let us assume that (26) holds in the exponential case, and let us consider the general case of a random variable X satisfying Assumption 1.1, then we have

E[p n,k (x)] = P(Λ(X) > Λ(a)|Λ(X) > Λ(x)) = exp(-Λ(a) + Λ(x)) = P(X > a|X > x) = P (x).
The first equality is a consequence of Corollary 3.4 and Theorem 4.1 in the exponential case. The third equality is a direct consequence of the definition (21) of Λ.

The aim of this section is thus to prove Theorem 4.1 in the exponential case. In all the following, we denote by f (x) = exp(-x) the density of X, and we will use the notation introduced in Section 3.1.2 above for the density of the k-th statistics. Actually, the proof given below is valid as soon as X has a density f : we do not use the specific form of the density, and this specific form would not make the argument easier.

The proof of this result is divided into two steps. First, we show that the function x → p n,k (x) is solution of a functional equation. Second, we show that the function x → P (x) is its unique solution.

4.1. Derivation of the functional equation satisfied by p n,k . Proposition 4.2. Let us assume that X admits a density f . Assume p = P (0) > 0. The function x ∈ [0, a] → p n,k (x) is solution of the following functional equation (with unknown q): for any 0 ≤ x ≤ a

(27) q(x) = a x 1 - k n q(y)f n,k (y; x)dy + θ n,k p (x), with (28) θ n,k p (x) = k-1 l=0 n -l n P S(x) n (l) ≤ a ≤ S(x) n (l+1) ,
where (S(x) n l ) 1≤l≤n are independent and identically distributed with density f (.; x) (see (23)), while for 1 ≤ l ≤ k -1, S(x) n (l) denotes the l-th order statistics of this n-sample. By convention, we set S(x) n (0) = x.

Proof:

The key idea is to decompose the expectation E[p n,k ] according to the (random but almost surely finite) value of the number J n,k (x) of iterations. The function θ n,k p appears as the result of the algorithm when J n,k (x) = 0, while the integral formulation corresponds to the case J n,k (x) > 0. In the latter case, we then condition on the value of the first level Z 1 = X 0 (k) and use Corollary 3.5. More precisely, we have

p n,k (x) = E pn,k (x) = E pn,k (x)1 J n,k (x)=0 + E pn,k (x)1 J n,k (x)>0 .
First, we have from ( 9) and ( 10)

E pn,k (x)1 J n,k (x)=0 = E C n,k (x)1 J n,k (x)=0 = E C n,k (x)1 X 0 (k) >a = k-1 l=0 n -l n E 1 X 0 (l) ≤a<X 0 (l+1) = θ n,k p (x).
Second, using conditional expectation with respect to Z 1 ,

E pn,k (x)1 J n,k (x)>0 = E E 1 - k n 1 - k n J n,k (x)-1 C n,k (x)|Z 1 1 Z 1 ≤a = E 1 - k n E 1 - k n J n,k (Z 1 ) C n,k (Z 1 )|Z 1 1 Z 1 ≤a = 1 - k n E E pn,k (Z 1 )|Z 1 1 Z 1 ≤a = 1 - k n E p n,k (Z 1 )1 Z 1 ≤a = 1 - k n a x p n,k (y)f n,k (y; x)dy,
where on the event

Z 1 ≤ a the equality E[(1 -k n ) J n,k (x)-1 C n,k (x)|Z 1 ] = E[(1 -k n ) J n,k (Z 1 ) C n,k (Z 1 )|Z 1
] is a consequence of Corollary 3.5, and of the fact that both J n,k (x) and J n,k (Z 1 ) are almost surely finite.

The intuition for this computation is that after the first step, if the algorithm does not stop, we just have to restart the algorithm from the level Z 1 , and consider its associated estimator of the probability. The multiplication by (1 -k n ) corresponds to the first iteration, which allows to go from level Z 0 = x to level Z 1 .

We will need in the following a more explicit formula for θ n,k p (x). Lemma 4.3. We have for any x ∈ [0, a]

(29) θ n,k p (x) = (1 -F (a; x)) (1 -F n-1,k (a; x)
) . Proof: We recall that (S(x) n l ) 1≤l≤n denotes a n-sample of i.i.d. random variables with law L(X|X > x), and that for any l ∈ {1, . . . , n} the random variable S(x) n (l) denotes the l-th order statistics of this n-sample: almost surely we have S(x) n (1) < . . . < S(x) n (n) ; by convention we moreover have S(x) n (0) = x. The proof is based on the partition of the n-sample into the (n -1)-sample (S(x) n-1 l ) 1≤l≤n-1 := (S(x) n l ) 1≤l≤n-1 and the random variable S(x) n n . We express the probabilities appearing in the definition of θ n,k p (x), using the cumulative distribution function of S(x) n n and of the l-th order statistics S(x) n-1

(l)
of the (n -1)-sample.

First, starting from (28), we write

k-1 l=0 n -l n P S(x) n (l) ≤ a ≤ S(x) n (l+1) = k-1 l=0 n -l n (n!) P S(x) n 1 ≤ . . . ≤ S(x) n l ≤ a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n = k-1 l=0 (n -l) ((n -1)!) P (S(x) n 1 ≤ . . . ≤ S(x) n l ≤ a) P a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n ,
since the random variables S(x) n l are independent, for l ∈ {1, . . . , n}. Now, for a fixed l ∈ {0, . . . , k -1}, using the fact that (S(x) n h ) l+1≤h≤n are i.i.d. and changing the position of S(x) n n in the ordered sample S(x) n l+1 ≤ . . . ≤ S(x) n n-1 , we have for any j ∈ {l, . . . , n -1}

P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n ) = P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n j ≤ S(x) n n ≤ S(x) n j+1 ≤ . . . ≤ S(x) n n-1 ),
with the convention that for j = l the right-hand side above is

P(a ≤ S(x) n n ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n-1 ), while for j = n -1 it is P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n-1 ≤ S(x) n n ).
We obtain (since all the terms in the sum below are all the same)

(n -l)P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n ) = n-1 j=l P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n j ≤ S(x) n n ≤ S(x) n j+1 ≤ . . . ≤ S(x) n n-1 ) = P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n-1 , S(x) n n ≥ a) = P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n-1 )P(S(x) n n ≥ a)
. The last equality comes from independence, and the sum expresses the fact that there are nl positions to insert S(x) n n in the increasing sequence S(x) n l+1 ≤ . . . ≤ S(x) n n-1 . Thus

k-1 l=0 n -l n P S(x) n (l) ≤ a ≤ S(x) n (l+1) = k-1 l=0 (n -1)! P S(x) n 1 ≤ . . . ≤ S(x) n l ≤ a)P(a ≤ S(x) n l+1 ≤ . . . ≤ S(x) n n-1 )P(S(x) n n ≥ a = P(S(x) n n ≥ a) k-1 l=0 (n -1)! P S(x) n-1 1 ≤ . . . ≤ S(x) n-1 l ≤ a ≤ S(x) n-1 l+1 ≤ . . . ≤ S(x) n-1 n-1 = P(S(x) n n ≥ a) k-1 l=0 P S(x) n-1 (l) ≤ a ≤ S(x) n-1 (l+1) = P(S(x) n n ≥ a)P(S(x) n-1 (k) ≥ a) = (1 -F (a; x)) (1 -F n-1,k (a; x)) .
This concludes the proof of Lemma 4.3.

Notice that we have proved a stronger statement: for l ∈ {0, . . . , k -1} Proof: Let p 1 , p 2 : [0, a] → R + be two bounded solutions. Then, we have for any x ∈ [0, a]

(30) n -l n P S(x) n (l) ≤ a ≤ S(x) n (l+1) = P(S(x) n n ≥ a)P S(x) n-1 (l) ≤ a ≤ S(x) n-1 (l+1) .
|p 1 (x) -p 2 (x)| ≤ 1 - k n a x |p 1 (y) -p 2 (y)|f n,k (y; x)dy ≤ 1 - k n p 1 -p 2 ∞ a x f n,k (y; x)dy ≤ 1 - k n p 1 -p 2 ∞ , which shows that p 1 = p 2 , since k/n > 0.
Notice that both functions p n,k and P take values in [0, 1], and are therefore bounded. Thanks to Proposition 4.2, p n,k satisfies (27), and Theorem 4.1 is thus a direct consequence of Lemma 4.4 if we prove that P is also solution of this functional equation.

Proof of Theorem 4.1:

The proof consists in proving that x → P (x) = 1 -F (a; x) is solution of (27). For this we have to compute for x ∈ [0, a],

a x 1 - k n P (y)f n,k (y; x)dy = a x (1 -F (a; y)) (n -k)k n n k F (y; x) k-1 f (y; x)(1 -F (y; x)) n-k dy = a x (1 -F (a; y)) (1 -F (y; x)) k n -1 k F (y; x) k-1 f (y; x)(1 -F (y; x)) n-k-1 dy = a x (1 -F (a; y)) (1 -F (y; x)) f n-1,k (y; x)dy = (1 -F (a; x)) a x f n-1,k (y; x)dy = (1 -F (a; x)) F n-1,k (a; x),
thanks to the definitions (23) of f n,k and f n-1,k , and the relation

(1 -F (a; y)) (1 -F (y; x)) = 1 -F (a; x)
for any x ≤ y ≤ a, obtained from (19).

We then conclude by checking the following identity, which is a consequence of Lemma 4.3.

a x 1 - k n P (y)f n,k (y; x)dy + θ n,k p (x) = (1 -F (a; x)) F n-1,k (a; x) + (1 -F (a; x)) (1 -F n-1,k (a; x)) = 1 -F (a; x) = P (x)
which concludes the proof.

Variance and computational cost: results

In this Section, we introduce a notion of cost for the algorithm (related to the variance of the estimator and to the expected number of iterations), which allows to study the influence of the parameters n and k. We then give asymptotic expansions of the variance, the expected number of iterations, and the cost, when n tends to +∞, for fixed values of k and of the probability p, and give interpretations of the results, compared to a direct Monte Carlo estimate as presented in the Introduction.

All these results are given under Assumption 1.1. The proofs are then given in Section 6.

5.1. Definition of the cost. In the following, we denote by c 0 the cost corresponding to the simulation of one random variable sampled according to the law L(X|X > x), for any x ∈ [0, a]. We assume that this cost does not depend on x.

We consider the Monte Carlo approximation of p using M independent realizations of the AMS algorithm 2.1. The associated estimator is

(31) p n,k M (x) = 1 M M m=1 pn,k m (x),
where pn,k m (x) is the AMS estimator for the m-th independent realization of the algorithm. Following the reasoning used in the introduction on the direct Monte Carlo estimator, for a given tolerance error ǫ, we want the relative error Var(p n,k M (x))

1/2 P (x) to be less than ǫ, i.e.

1 M Var(p n,k (x)) ≤ ǫ 2 P (x) 2 .
We thus have to choose M = Var( pn,k (x)) ǫ 2 P (x) 2 . For each realization m of the algorithm, J n,k m (x) iterations are necessary, so that k J n,k m (x) + n random variables are sampled: n at the initial step, and then k new ones at each iteration. This gives a cost c 0 k J n,k m (x) + n , where c 0 is the computational cost of the sampling of one random variable distributed according to L(X|X > x). Moreover, at the first iteration, one needs to sort the random variables (X 0 1 , . . . , X 0 n ) (with an associated cost c 1 n log n) and at each iteration, one has to insert the k new sampled particles into the already sorted (nk) remaining particles (with an associated cost c 1 k log n). The sorting procedures are thus associated with a cost c 1 (log n) k J n,k m (x) + n . The total cost is thus

M m=1 (c 0 + c 1 log n) k J n,k m (x) + n
and by an application of the Law of Large Numbers, when M is large, it is legitimate to consider that the cost to obtain a relative error of size ǫ is thus

(32) (c 0 + c 1 log n) C n,k (x) ǫ 2 where C n,k (x) = Var(p n,k (x)) k E[J n,k (x)] + n P (x) 2 .
This is consistent with the standard definition of the efficiency of a Monte Carlo procedure as "inversely proportional to the product of the sampling variance and the amount of labour expended in obtaining this estimate", see [START_REF] Hammersley | Monte Carlo methods[END_REF]Section 2.5]. This should be compared with the cost of a direct Monte Carlo computation, which we recall (see (4)):

(33)

(1 -P (x))

ǫ 2 P (x) c 0 .
Remark 5.1. If we add the possibility of using N ≥ 1 processors to sample in parallel the required random variables at each iteration (assuming for simplicity that k/N is an integer) then the cost is divided by N , and is thus

(c 0 + c 1 log n) C n,k (x) ǫ 2 N .
Notice that this resulting cost is the same as if we run in parallel N independent realizations of pn,k m (x) to compute the estimator p n,k M (x) (assuming for simplicity that M/N is an integer). Both these parallelization strategies have the same effect on the cost in the setting of this article.

Let us set a few notations: for x ∈ [0, a],

(34) v n,k (x) = E (p n,k (x)) 2 and T n,k (x) = E J n,k (x) + 1.
Notice that T n,k (x) is the expected number of steps in the algorithm (the initialization plus J n,k (x) iterations). Using this notation, we have

(35) C n,k (x) = v n,k (x) -P (x) 2 P (x) 2 k T n,k (x) + n -k .
5.2. Asymptotic expansions of the variance and of the computational cost. We divide the results into three parts. We first study the variance, and then the average number of iterations. Finally, we combine the results to get the cost. We do not have explicit expressions for each value of k and n, but we get informative asymptotic results when n → +∞. We assume that p = P (0) > 0, and consider x ∈ [0, a), such that P (x) < 1. Note that if P (x) = 1, then almost surely p n,k (x) = 1, Var(p n,k (a)) = 0 and T n,k (a) = 1, so that no asymptotic expansion is necessary. Proposition 5.2. For any fixed k and any 0 ≤ x < a with P (x) < 1, when n → +∞ we have

(36) Var(p n,k (x)) = P (x) 2 n   -log(P (x)) + ((log(P (x))) 2 -log(P (x)) (k -1) 2n + o 1 n   .
Proposition 5.3. For any fixed k and 0 ≤ x < a with P (x) < 1, when n → +∞ we have

(37) T n,k (x) = n -log(P (x)) 1 k - k -1 2kn + 3k -1 2kn + o 1 n .
Finally, we have the following result on the cost:

Theorem 5.4. For any fixed k and 0 ≤ x < a with P (x) < 1, when n → +∞ we have

(38) C n,k (x) = (log(P (x))) 2 -log(P (x)) + 1 n -log(P (x)) [k -1] + 1 2 (log(P (x))) 2 - 1 2 (log(P (x))) 3 + o 1 n .
The proof of Theorem 5.4 from the Propositions 5.2 and 5.3 is straightforward using (35). The proof of the two Propositions is long and technical, and is postponed to Section 6.

Let us also state an immediate corollary of Theorem 4.1 and of Proposition 5.2.

Corollary 5.5. For any fixed k and 0 ≤ x < a with P (x) < 1, when n → +∞, we have the following convergence in probability:

pn,k (x) → P (x).

Remark 5.6. The case k = 1 is much simpler than the general case k > 1, and allows for direct computations. As seen in Section 3.3, for n ≥ 2 and 0 ≤ x ≤ a, J n,1 (x) has a Poisson distribution with parameter -n log(P (x)). We can then easily check the unbiased property p n,1 (x) = P (x), and compute T n,1 (x) = -n log(P (x)) + 1 and Var(p n,k (x)) = P (x) 2 (P (x) -1/n -1). In particular, no asymptotic expansions are required in order to understand the behavior of the computational cost. Notice that in the case k = 1, E J n,1 (x) = -n log(P (x)), and thus J n,1 (

n is an unbiased estimator of log(P (x)), with variance -log(P (x)) n . There is no such statement when k > 2: from Proposition 5.3, we have the limit

E J n,k n → -log(P (x))
k , but the first-order term in 1/n is equal to k-1 2k (1 + log(P (x))) and is therefore not zero, except if k = 1 or P (x) = exp(-1) -in which case one should compute an higher-order expansion to prove that there is a bias.

Let us now recall the main two practical consequences of Theorem 5.4, that we already stressed in the introduction. First, in view of (33) and (32)-(38), the AMS algorithm is more efficient than a direct Monte Carlo procedure to estimate p = P (0) if

1 -p p c 0 > (c 0 + c 1 log(n)) (log p) 2 -log p
which is always true for sufficiently small p. Second, from Theorem 5.4, we observe that all choices of k give the same leading order term for the cost. But looking at the term of order 1 n , we see that the optimal choice is k = 1. This conclusion can also be deduced from the asymptotic expansion on the variance given in Proposition 5.2.

Proof of the variance and computational cost estimates

This Section is devoted to the proof of Propositions 5.2 and 5.3, which together yield the cost estimate of Theorem 5.4. The main steps are the following:

• In Section 6.1, we first show that we can reduce the analysis to the exponential case, by the change of variable using the function Λ, as explained in Section 3.2. • In Section 6.2, we then derive functional equations on the two functions v n,k and T n,k (defined by (34)). These functional equations can actually be obtained not only in the exponential case, but for any X which admits a density with respect to the Lebesgue measure. • In Section 6.3, we prove that, in the exponential case, the functional equations on v n,k and T n,k are equivalent to linear ordinary differential equations of order k. • Finally, we compute asymptotic expansions of v n,k (in Section 6.4) and T n,k (in Section 6.5) in the large n limit, for fixed k and p.

The main simplification provided by the exponential case is that the functional equations can be recast as linear differential equations (see Remark 6.5 below).

6.1. Reduction to the exponential case. We have seen in Section 3.2 (see Corollary 3.4) that the estimator pn,k (x) obtained with the AMS algorithm applied to a general random variable X (satisfying Assumption 1.1) with initial condition x and target level a is exactly the same in law as the estimator pn,k (Λ(x)) that is obtained with the AMS algorithm applied to an exponentially distributed random variable X, with initial condition Λ(x) and target level Λ(a).

It is therefore sufficient to prove the Propositions 5.2 and 5.3 in the exponential case. Indeed, if we obtain in the exponential case, for an initial condition x and a target level a ≥ x:

Var(p n,k (x)) = exp(2(x -a)) n   (a -x) + (a -x) 2 + (a -x) (k -1) 2n + o 1 n   , (39)
T n,k (x) = n (a -x) 1 k - k -1 2kn + 3k -1 2kn + o 1 n , (40)
then the general case is easily obtained by replacing x by Λ(x) and a by Λ(a), since Λ(x)-Λ(a) = -Λ(a; x) = log(1 -F (a; x)) = log(P (x)).

6.2. Functional equations satisfied by v n,k and T n,k . We now write functional equations satisfied by v n,k (x) and T n,k (x) defined by (34). Even though we will only need these functional equations in the exponential case as explained above, we derive these functional equations in a more general setting, namely when X admits a density f with respect to the Lebesgue measure. We refer to Section 3.1.2 for relevant notation associated to this setting. Notice that the derivations of these functional equations are very similar to the derivation of the functional equation on p n,k in the proof of Proposition 4.2. Proposition 6.1. Assume P (0) > 0. The function x → v n,k (x) is solution of the following functional equation (with unknown w):

for any 0 ≤ x ≤ a (41) w(x) = a x 1 - k n w(y)f n,k (y; x)dy + θ n,k v (x) with (42) θ n,k v (x) = k-1 l=0 (n -l) 2 n 2 P S(x) n (l) ≤ a ≤ S(x) n (l+1) .
where (S(x) n l ) 1≤l≤n are independent and identically distributed with density f (x, .), and S(x) n (l) denotes the l-th order statistics of this n-sample. By convention, S(x) n (0) = x.

Similarly, we can derive a functional equation satisfied by T n,k .

following Linear Ordinary Differential Equations (ODEs) of order k: for x ∈ [0, a]:

d k dx k p n,k (x) = 1 - k n µ n,k p n,k (x) + k-1 m=0 r n,k m d m dx m p n,k (x), (45)
d k dx k v n,k (x) = 1 - k n 2 µ n,k v n,k (x) + k-1 m=0 r n,k m d m dx m v n,k (x), (46)
d k dx k T n,k (x) = k-1 m=1 r n,k m d m dx m T n,k (x) + µ n,k . (47)
Notice that in (47) the summation starts at m = 1, while in (45) and (46) it starts at m = 0. The coefficients µ k,n and (r k,n m ) 0≤m≤k-1 are defined by a simple induction formula, see (55). Moreover, the functions p n,k , v n,k and T n,k satisfy the following boundary conditions at point x = a: for m ∈ {0, . . . , k -1}

d m dx m p n,k (x) x=a = 1, (48)
d m dx m v n,k (x) x=a = 1 n + 1 - 1 n 2 m , (49)
d m dx m T n,k (x) x=a = 1 m=0 . ( 50
)
The main tool for the proof of Proposition 6.4 is the following formula on the derivative of the density f n,k (.; x) with respect to x: for all y > x, (51)

d dx f n,1 (y; x) = nf n,1 (y; x) for k ∈ {2, . . . , n -1}, d dx f n,k (y; x) = (n -k + 1)(f n,k (y; x) -f n,k-1 (y; x)).
Recall that f (y; x) = exp(-(yx)) for y ≥ x. The proof of the first formula in (51) is straightforward, since f n,1 (y; x) = n exp(-n(yx)). For k ∈ {2, . . . , n -1}, we write (using (23))

d dx f n,k (y; x) = d dx k n k F (y; x) k-1 f (y; x)(1 -F (y; x)) n-k = k n k d dx (1 -exp(x -y)) k-1 exp ((n -k + 1)(x -y)) = k n k -(k -1) exp(x -y)(1 -exp(x -y)) k-2 exp ((n -k + 1)(x -y)) + (n -k + 1)(1 -exp(x -y)) k-1 exp ((n -k + 1)(x -y)) = (n -k + 1)f n,k (y; x) -(k -1) k n k (k -1) n k-1 f n,k-1 (y; x) = (n -k + 1) (f n,k (y; x) -f n,k-1 (y; x)) .
Remark 6.5. A generalization of (51) holds in a more general case than the exponential setting. Indeed, if X has a density f , then, for y > x (52)

d dx f n,1 (y; x) = f (x) 1 -F (x)
nf n,1 (y; x)

for k ∈ {2, . . . , n -1}, d dx f n,k (y; x) = f (x) 1 -F (x) (n -k + 1) (f n,k (y; x) -f n,k-1 (y; x)) .
In the exponential case, the simplification f (x) 1-F (x) = 1 helps getting simpler formulae, which lead to the linear ODEs of Proposition 6.4. It is also worth noting the following formula f (x) 1-F (x) = -d dx log(1-F (x)) = d dx Λ(x), which explains the role played by the change of variable using the function Λ to reduce the general case to the exponential case.

Proof of Proposition 6.4: We mainly focus on the derivation of the ODE (45) for p n,k . The ODEs (46) and (47) are obtained with similar arguments.

For any 1 ≤ l ≤ k, we define for 0 ≤ x ≤ a (53)

I n,k l (x) = a x 1 - k n p n,k (
y)f n,l (y; x) dy.

We also set

I n,k 0 (x) = (1 -k n )p n,k (x). As a consequence of (51), we get for 1 ≤ l ≤ k d dx I n,k l (x) = (n -l + 1)(I n,k l (x) -I n,k l-1 (x)).
Precisely, for 2 ≤ l ≤ k, this formula directly follows from (53) and f n,l (x; x) = 0. When l = 1,

d dx I n,k 1 (x) = -1 - k n p n,k (x)f n,1 (x; x) + nI n,k 1 (x) = n(I n,k 1 (x) -I n,k 0 (x)).
The ODE (45) on p n,k is then obtained as follows.

• p n,k (x)θ n,k p (x) = I n,k k (x) (this is the functional equation ( 27)). • For any 0 ≤ l ≤ k, we prove by induction that the following formula holds:

(54) d l dx l p n,k (x) -θ n,k p (x) = µ n,k l I n,k l (x) + l-1 m=0 r n,k m,l d m dx m p n,k (x) -θ n,k p (x) .
The coefficients are defined recursively as follows:

(55)

µ n,k 0 = 1, µ n,k l+1 = -(n -k + l + 1)µ n,k l ;      r n,k 0,l+1 = -(n -k + l + 1)r n,k 0,l , if l > 0, r n,k m,l+1 = r n,k m-1,l -(n -k + l + 1)r n,k m,l , 1 ≤ m ≤ l -1, r n,k l,l+1 = n -k + l + 1.
• The ODE is obtained at l = k, using the definition of I n,k 0 and setting µ n,k = µ n,k k and r n,k m = r n,k m,k :

d k dx k p n,k (x) -θ n,k p (x) = 1 - k n µ n,k p n,k (x) + k-1 m=0 r n,k m d m dx m p n,k (x) -θ n,k p (x) .
Similarly, we obtain

d k dx k v n,k (x) -θ n,k v (x) = 1 - k n 2 µ n,k v n,k (x) + k-1 m=0 r n,k m d m dx m v n,k (x) -θ n,k v (x) and
d k dx k T n,k (x) = µ n,k T n,k (x) -r n,k 0 + k-1 m=0 r n,k m d m dx m T n,k (x). Since r n,k 0 = (-1) k-1 n . . . (n -k + 1) = -µ n,k
, after simplifications we obtain (47).

To obtain the ODEs given in Proposition 6.4, it remains to prove that (56)

d k dx k θ n,k p (x) = k-1 m=0 r n,k m d m dx m θ n,k p (x), d k dx k θ n,k v (x) = k-1 m=0 r n,k m d m dx m θ n,k v (x).
The argument is as follows. Using Lemma 4.3, in the special case of exponential random variables, elementary computations show that for 0

≤ x ≤ a θ n,k p (x) = k-1 j=0 k n -1 k k -1 j (-1) j n -k + j exp ((n -k + j + 1)(x -a)) .
This shows that θ n,k p as well as its derivatives, are linear combinations of the linearly independent functions x → exp((nk + 1)x), . . ., x → exp(nx). For our purpose, the exact expression of the coefficients does not matter.

But from Theorem 4.1, we know that p n,k (x) = P (x) = exp(xa). We thus conclude by a linear independence argument that

d k dx k p n,k (x) = 1 - k n µ n,k p n,k (x) + k-1 m=0 r n,k m d m dx m p n,k (x), d k dx k θ n,k p (x) = k-1 m=0 r n,k m d m dx m θ n,k p (x).
In particular, the second equality implies that for 0 ≤ j ≤ k -1

(57) d k dx k exp ((n -k + j + 1)(x -a)) = k-1 m=0 r n,k m d m dx m exp ((n -k + j + 1)(x -a)) .
The second equality of (56) can be proven using the same arguments.

Remark 6.6. It is actually possible to prove directly the identity (57) from the definition of the coefficients r n,k m , without resorting to the result of Theorem 4.1. Indeed, let us introduce S j l := (nk + j + 1) ll-1 m=0 r n,k m,l (n-k+j +1) m . Thanks to the recursion formula (55), introducing an appropriate telescoping sum, one obtains that the coefficients (S j l ) satisfy S j 0 = 1 and S j l+1 = (jl)S j l which implies that S j k = S j j+1 = 0 for 0 ≤ j ≤ k -1. One thus obtains that p n,k satisfies the ODE (45) without using Theorem 4.1.

This approach actually yields an alternative proof of Theorem 4.1 in the exponential case, that can then be extended to the general case using the reduction to the exponential case explained in Section 3.2. For this, we check that x → exp(xa) is solution of the differential equation on p n,k . First, it satisfies the appropriate condition at x = a given in Proposition 6.4. Second, we observe that t = 1 is a root of the characteristic polynomial equation associated with the linear ODE which writes (the unknown variable being t):

(n -t) . . . (n -k + 1 -t) n . . . (n -k + 1) = (n -k) n .
See (59) below for the derivation of the caracteristic equation.

To conclude the proof of Proposition 6.4, it remains to show the boundary conditions at x = a. From the recursion equation (54) leading to the differential equation on p n,k , and Lemma 4.3, we have for 0

≤ l ≤ k -1 d l dx l p n,k (x) x=a = d l dx l θ n,k p (x) x=a = d l dx l (exp(x -a)(1 -F n-1,k (a; x))) x=a = 1 + l j=1 l j d j dx j (1 -F n-1,k (a; x)) x=a = 1 - l j=1 l j d j-1 dx j-1 f n-1,k (a; x) x=a = 1.
Indeed, the second identity of (51) yields

d j-1 dx j-1 f n-1,k (a; x) x=a = 0 for any 1 ≤ j ≤ k -1.
Similarly, from the derivation of the ODE on v n,k we have

d l dx l v n,k (x) x=a = d l dx l θ n,k v (x) x=a = 1 n d l dx l (exp(x -a) (1 -F n-1,k (a; x))) x=a + 1 - 1 n d l dx l (exp(2(x -a)) (1 -F n-2,k (a; x))) x=a = 1 n + 1 - 1 n 2 l ,
using Lemma 6.3, and the same arguments as for p n,k . Finally, using similar arguments,

T n,k (a) = 1;

d l dx l T n,k (x) x=a = d l dx l 1 x=a = 0, 1 ≤ l ≤ k -1.
This concludes the proof of Proposition 6.4.

In the next Sections, we analyze the differential equations on v n,k and T n,k . We are not able to derive explicit expressions for the solutions, except when k = 1 (see Remark 5.6). However, we are able to analyze quantitatively the behavior when n → +∞.

6.4. Asymptotic expansion for the variance. In this Section, we prove Proposition 5.2 in the exponential case, namely (39). Since

Var(p n,k (x)) = v n,k (x) -(p n,k (x)) 2 ,
and we know from Theorem 4.1 that p n,k (x) = P (x) = exp(xa), we focus on v n,k . This function is solution of the linear ODE of order k given in Proposition 6.4, which we rewrite here:

d k dx k v n,k (x) = 1 - k n 2 µ n,k v n,k (x) + k-1 m=0 r n,k m d m dx m v n,k (x).
To understand how the solution v n,k behaves, we thus need to study the following associated polynomial equation with unknown t:

t k - k-1 m=0 r n,k m t m -µ n,k 1 - k n 2 = 0.
to (49)). Using Cramer's rule, we give explicit formulae and then get asymptotic expansions for each η l n,k when n → +∞. Proposition 6.8. Let k be fixed. When n → +∞, we have

η 1 n,k = 1 + o 1 n η l n,k = o(1)
, for k ∈ {2, . . . , k}. Thanks to the previous result and to the expression of v n,k as a combination of exponential functions (60), subtracting P (x) 2 = exp (2(xa)), we obtain the desired asymptotic expansion (39) for the variance when n → +∞. Notice that since x is assumed to be strictly smaller than a, the only contribution which remains is η 1 n,k exp β 1 n,k (xa) . The other terms in the sum (60) vanish exponentially fast. This concludes the proof of Proposition 5.2.

We end this Section with the proof of Proposition 6.8. Proof of Proposition 6.8: The family (η l n,k ) 1≤l≤k is solution of the following system of linear equations (using (49)):

             η 1 n,k + η 2 n,k + . . . + η k n,k = v n,k (a) = 1 η 1 n,k β 1 n,k + η 2 n,k β 2 n,k + . . . + η k n,k β k n,k = d dx v n,k (x) x=a = 2 -1 n • • • η 1 n,k (β 1 n,k ) k-1 + η 2 n,k (β 2 n,k ) k-1 + . . . + η k n,k (β k n,k ) k-1 = d k-1 dx k-1 v n,k (x) x=a = 1 n + 1 -1 n 2 k-1
. Using Cramer's rule (which gives the solution of an invertible linear system thanks to ratios of determinants), we see that

η 1 n,k = 1 n V (1, . . . , β k n,k ) V (β 1 n,k , . . . , β k n,k ) + 1 - 1 n V (2, . . . , β k n,k ) V (β 1 n,k , . . . , β k n,k )
, where V (λ 1 , . . . , λ k ) = det(λ i-1 j ) 1≤i,j≤k denotes the Vandermonde determinant of the complex numbers (λ i ) 1≤i≤k . We recall that V (λ 1 , . . . , λ k ) = 1≤i<j≤k (λ jλ i ).

Straightforward simplifications then imply that .

The asymptotic results on the coefficients β l n,k given in Proposition 6.7 finally show that η 1 n,k = 1 + o 1 n . The proof that η l n,k = o(1) for k ∈ {2, . . . , k} follows the same lines. 6.5. Asymptotic expansion for T n;k . In this Section, we prove Proposition 5.3 in the exponential case, namely (40), following the same approach as in the previous Section. Recall from Proposition 6.4 that T n,k is solution of the following linear differential equation of order k: This formulation allows to get the following analog of Proposition 6.7, on the k roots (α l n,k ) 1≤l≤k ∈ C k of (61). Proposition 6.9. Let k be fixed. When n → +∞, the roots (α l n,k ) 1≤l≤k of (61) satisfy: α 1 n,k = 0 α l n,k ∼ n (1exp(i2π(l -1)/k)) , for l ∈ {2, . . . k}.

d k dx k T n,k (x) -
We omit the proof of Proposition 6.9, since it is very similar to the proof of Proposition 6.7. We have already identified that α 1 n,k = 0 and the asymptotic formulae for α l n,k when l ∈ {2, . . . , k} are obtained with exactly the same arguments as for β l n,k in Proposition 6.7. Again, for n large enough, the roots (α l n,k ) l∈{1,...,k} are pairwise distinct.

The two differences with the analysis performed in the previous Section are the following. First, the differential equation (47) on T n,k contains a non-zero right-hand side (namely a constant). Moreover, constant functions are solutions of the differential equation without this right-hand side, since 0 is a root of (61).

Therefore T n,k can be expressed as a sum of an affine function and of exponential functions, for n large enough (so that roots are pairwise distinct): for any 0 ≤ x ≤ a (62)

T n,k (x) = ∆ n,k (ax) + δ 1 n,k + k l=2 δ l n,k exp(α l n,k (xa)),

for some complex coefficients ∆ n,k and δ l n,k , for 1 ≤ l ≤ k. We prove the following result, which then yields Proposition 5.3: Proposition 6.10. Let k be fixed. When n → +∞, we have zf n,k (z; 0)dz. Then, using (51) and an integration by parts, we get R( 1) . If we now write that z k -1 = R(z)(z -1), differentiating and setting z = 1 we obtain R(1) = k; a similar argument shows that R ′ (1) = k(k-1)

∆ n,k = n 1 k - k -1 2kn + o 1 n , δ 1
M n,1 = 1 n , M n,k = M n,k-1 + 1 n -k + 1 .

2

, so that finally

δ 1 ∞,k = 1 + k -1 2k = 3k -1 2k .
A similar Cramer's rule holds for each δ l n,k when l ≥ 2. It is then easy to check that δ l n,k = O(1). In fact, an analytical formula for the limit of δ l n,k in the limit n → +∞ can be written, but we do not need such a sharp result for our purposes.

This concludes the proof of Proposition 6.9.

Definition 2 . 3 .

 23 Let A and B be two disjoint closed subsets in R d . The committor function ξ associated with the SDE (12) and the sets A and B is the unique solution of the following partial differential equation (PDE):

2 . 1 .Lemma 3 . 1 .

 2131 The basic remark is the following. Since by Assumption 1.1 the cumulative distribution function F (.; x) is continuous, we have the following classical result, see for instance Proposition 2.2 in[START_REF] Asmussen | Stochastic simulation: Algorithms and analysis Springer[END_REF]. If Y ∼ L(X|X > x), then F (Y ; x) is uniformly distributed on (0, 1), and thus Λ(Y ; x) has an exponential law with parameter 1.

A = P Y j+1 1 -Y j+1 1 -

 11 . n} = I -⊔ I + , such that Card I -= k. S j+1 > y 1 , . . . , Y j+1 n -S j+1 > y n , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 = I -⊂{1,...,n} Card(I-)=k P S j+1 > y 1 , . . . , Y j+1 n -S j+1 > y n , S 1 -S 0 > s 1 , . . . , S j+1 -S j > s j+1 , σ j+1 (1), . . . , σ j+1 (k) = I - = I -⊂{1,...,n} Card(I-)=k

4. 2 .Lemma 4 . 4 .

 244 Proof of Theorem 4.1. Let us first state a uniqueness result. The functional equation (27) admits at most one solution p : [0, a] → R + in L ∞ ([0, a]).

  m T n,k (x) = µ n,k .The associated polynomial equation ykk-1 m=1 r n,km y m = 0 admits 0 as a root. Moreover, thanks to (58) it can be rewritten as(61) (nt) . . . (nk + 1t) n . . . (nk + 1) = 1.

  = O(1), l ∈ {2, . . . k}. Proof: First, the differential equation (47) on T n,k in Proposition 6.4 is in fact valid on (-∞, a], not only on [0, a]. It corresponds to the estimation by the AMS algorithm of P (x) = P(x + X > a) = exp(xa) for x ≤ a and X ∼ E(1) exponentially distributed.Inserting (62) into the functional equation (43) and letting x → -∞ yields∆ n,k = 1 +∞ 0 z f n,k (z; 0)dz .Indeed, we get, using the change of variable z = yx and the fact that f n,k (y; x) = f n,k (yx; 0)T n,k (x) = 1 + a x ∆ n,k (xy)f n,k (y; x)dy + a x ∆ n,k (ax)f n,k (y; x)dy n,k (y-a) f n,k (y; x)dy = 1 -∆ n,k a-x 0 z f n,k (z; 0)dz + ∆ n,k (ax) z f n,k (z; 0)dz.Taking the limit x → -∞ above and in (62), we obtainT n,k (x) = ∆ n,k (ax) + δ 1 n,k + o(1), T n,k (x) = ∆ n,k (ax) + δ 1 n,k + 1 -∆ n,k a-x 0 z f n,k(z; 0)dz + o(1), and thus 1 -∆ n,k +∞ 0 f n,k (z; 0)dz = 0. Let us define M n,k = +∞ 0

( 1 - 1 kR ′ ( 1 )

 111 z) is a Vandermonde determinant, we have the explicit formulaQ(z) = V (-z, 1ξ 2,k , . . . , 1ξ k,k ) = V (1ξ 2,k , . . . , 1ξ k,k ) k l=2 ξ l,k + z) = V (1ξ 2,k , . . . , 1ξ k,k )R(z + 1), with R(z) = k l=2 (zξ l,k ) = z k -1 z-1 ; thus δ 1 ∞,k = 1 +

A precise mathematical statement is as follows. Let (U j ℓ ) 1≤ℓ≤k,j∈N * be i.i.d. random variables, uniformly distributed on (0, 1) and independent from all the other random variables. Then set χ j ℓ = F (.; Z j ) -1 (U j ℓ ), where F (.; x) -1 is the inverse distribution function associated with the (conditional) probability distribution L(X|X > x), see 19. We slightly abuse notation by using L(X|X > Z j ) rather than L(X|X > z)| z=Z j .

The uniqueness of the permutation σ j+1 is justified by Proposition

3.2.

Proposition 6.2. Assume P (0) > 0. The function x → T n,k (x) is solution of the following functional equation (with unknown T ): for any 0 ≤ x ≤ a (43) T (x) = a x T (y)f n,k (y; x)dy + 1.

We do not give the details of the proofs of these two results. The proof of Proposition 6.1 follows exactly the same lines as the proof of Proposition 4.2. For Proposition 6.2, using the same arguments, one obtains: Similarly to Lemma 4.3 on θ n,k p , we have a more explicit formula for θ n,k v that will be useful below.

Proof of Lemma 6.3: The notations are the following: S(x) n 1 , . . . , S(x) n n are n independent random variables distributed according to L(X|X > x), S(x) n-1

is the l-th order statistic of the random variables S(x) n-1

is the l-th order statistic of the random variables

Here again S(x) n (0) = S(x) n-1 (0) = x. From (42), and using twice the equality (30) obtained in the proof of Lemma 4.3, we have

As mentioned above, the functional equations ( 41)-( 43) and the equation (44) on θ n,k v are valid for any X with a density f . However, we are only able to exploit them in the exponential case. From now on, we thus only consider the exponential case:

6.3. Ordinary differential equations on p n,k , v n,k and T n,k in the exponential case. From the functional equations ( 27), ( 41) and (43), we show that the functions p n,k , v n,k and T n,k on [0, a] are solutions of linear ordinary differential equations, in the exponential case. Proposition 6.4. Let n and k ∈ {1, . . . , n -2} be fixed and let us assume that X ∼ E(1). There exist real numbers µ k,n and (r k,n m ) 0≤m≤k-1 , depending only on n and k, such that p n,k , v n,k and T n,k satisfy the

In order to study the behavior of the (complex) roots of this equation, we first observe that (58)

since by (57) the k roots of this polynomial function are n, . . ., nk + 1. Moreover, µ n,k = (-1) k n . . . (nk + 1); therefore, the polynomial equation can be rewritten

Proposition 6.7. Let k be fixed. There is a unique root of (59) in the real interval [START_REF] Asmussen | Stochastic simulation: Algorithms and analysis Springer[END_REF][START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF], denoted by β 1 n,k . The other roots of (59) in C are denoted by β 2 n,k , . . . , β k n,k .Moreover, we have the following asymptotic expansions when n → +∞:

Proof of Proposition 6.7: We observe that the complex numbers

are solutions of the following polynomial equation of degree k (with unknown t):

The claim for l ∈ {2, . . . , k} then follows by continuity of the roots of a polynomial function of fixed degree with respect to the coefficients. Indeed, in the limit n → +∞ we obtain the equation

. This argument for l = 1 yields that β 1 n,k /n goes to 0, and this is not sufficient for our purposes. To study the behavior of β 1 n,k , let us introduce the polynomial function

This function is strictly non-decreasing on the interval (-∞, nk + 1] (which contains no root of P n,k and of its derivative). Now straightforward computations show that P n,k (1)

, and thus P n,k admits a root in the interval [START_REF] Asmussen | Stochastic simulation: Algorithms and analysis Springer[END_REF][START_REF] Bréhier | Mathematical analysis of adaptive multilevel splitting algorithms[END_REF]. For n sufficiently large, since

Moreover, elementary computations show that

and therefore since P n+1,k is non-decreasing, β 1 n,k < β 1 n+1,k : the sequence (β n,k ) n>k is non-decreasing. Since β n,k ≤ 2 for any n > k, the sequence converges.

To identify the limit, and higher order terms in the expansion of β 1 n,k , one postulates an ansatz

For n large enough, all the roots are therefore simple, and we can express the function v n,k in the following way:

for some complex numbers (η l n,k ) 1≤l≤k , satisfying appropriate conditions to satisfy the boundary conditions (49). In particular, these complex numbers are such that v n,k is real-valued which implies necessarily η 1 n,k ∈ R. Actually, these complex numbers are solution to a system of linear equations (which corresponds Indeed,

We therefore obtain

The values of δ l n,k for 1 ≤ l ≤ k depend only on the conditions at x = a for T n,k and its derivatives up to order k -1, namely (50). They are solutions of a system of linear equations and they can be expressed thanks to Cramer's rule. The family (δ l n,k ) 1≤l≤k is solution of the following system of linear equations:

Using Cramer's rule and Vandermonde determinants, we see that

Let us introduce a few more notations. We define ξ l,k = exp(i2π(l -1)/k) for 2 ≤ l ≤ k, and the polynomial function:

(1ξ 2,k ) 2 . . .

Then by considering the limit n → +∞ in (63), plugging the asymptotic expansions for α l n,k and for ∆ n,k from Propositions 6.9 and 6.10, one obtains