Central Limit Theorem for Adaptative Multilevel Splitting Estimators in an Idealized Setting

0.75 MB·English

by  Charles-Edouard Bréhier

Central Limit Theorem for Adaptive Multilevel Splitting
Estimators in an Idealized Setting
Charles-Edouard Br´ehier ∗ Ludovic Gouden`ege † Lo¨ıc Tudela ‡
January 8, 2015
5
Abstract
1
0 TheAdaptiveMultilevelSplittingalgorithmisaverypowerfulandversatileiterativemethod
2 to estimate the probability of rare events, based on an interacting particle systems. In [3], in a
so-called idealized setting, the authors prove that some associated estimators are unbiased, for
n
a each value of the size n of the systems of replicas and of a resampling number k.
J Here we go beyond and prove these estimator’s asymptotic normality when n goes to +∞, for
7 anyfixedvalueofk. Themainingredientistheasymptoticanalysisofafunctionalequationon
anappropriatecharacteristicfunction.Somenumericalsimulationsillustratetheconvergenceto

] rely on Gaussian confidence intervals.
R
Key words: Monte-Carlosimulation,rareevents,multilevelsplitting,centrallimittheorem
P Subject Classifications: 65C05; 65C35; 60F05
.
h
at 1 Introduction
m
Manymodelsfromphysics,chemistryorbiologyinvolvestochasticsystemsfordifferentpurposes:
[
taking into account some uncertainty with respect to some data parameters, or to allow for

1 dynamical phase transitions between different configurations of the system. This phenomenon
v
often referred to as metastability is observed for instance when one uses a d-dimensional
9
overdamped Langevin dynamics
9
3 (cid:112)
dX =−∇V(X )dt+ 2β−1dW
t t t
1
0 associated with a potential function V with several local minima. Here W denotes a d-
.
1 dimensionalstandardWienerprocess. Whentheinversetemperatureβincreases,thetransitions
0 become rare events (their probability decreases exponentially fast).
5 Inthispaper,weadoptanumericalpointofview,andanalyzeamethodwhichoutperforms
1 a pure Monte-Carlo method for a given computational effort in the small probability regime
v: (in terms of relative error). Two important families of methods have been introduced in the
i 1950s and extensively developed later in order to efficiently address this rare event estimation
X
problem: importance sampling, and importance/multilevel splitting - see [11] and [9] for a
r more recent treatment. For more general references, we refer for instance to [12].
a
The method we study in this work belongs is a multilevel splitting algorithms. The main
advantage of this kind of methods is that they are non-intrusive: one does not need to modify
the model to go beyond a pure Monte-Carlo approximation scheme. The present method has
an additional feature: the so-called levels are computed adaptively. To explain more precisely
the algorithm and its properties, from now on we only focus on a simpler setting, though
paradigmatic of the rare event estimation problem.
Let X be a real random variable, and a some fixed, given threshold. We want to compute
an approximation of the tail probability p:=P(X >a). The splitting strategy in the regime
∗Universit´eParis-Est,CERMICS(ENPC),6-8-10AvenueBlaisePascal,Cit´eDescartes,F-77455Marne-la-Vall´ee,
France. INRIAParis-Rocquencourt,DomainedeVoluceau-Rocquencourt,B.P.105-78153LeChesnay,France.
e-mail:[email protected],e-mail:[email protected]
†F´ed´erationdeMath´ematiquesdel’E´coleCentraleParis,CNRS,Grandevoiedesvignes,92295Chaˆtenay-Malabry,
France. e-mail:[email protected]
‡EnsaeParisTech,3AvenuePierreLarousse,92240Malakoff,France. e-mail:[email protected]
1
when a becomes large consists in introducing the following decomposition of p as a product of
conditional probabilities:
P(X >a)=P(X >a |X >a )...P(X >a |X >a )P(X >a ),
n n−1 2 1 1
for some sequence of levels a < ... < a < a = a. The common interpretation of this
1 n−1 n
formula is that the event that X >a is split in n (conditional) crossing probabilities for X,
each much higher than p and thus easier to approximate (each independently).
Tooptimizethevariance,thelevelsmustbechosensuchthatalltheconditionalprobabilities
areequaltop1/n,withnaslargeaspossible. However,theapriori knowledgeoflevelssatisfying
this condition is not available in practical cases.
Notice that in principle to apply this splitting strategy, one needs to know how to sample
according to the conditional distributions appearing in the splitting formula. If this condition
holds, we say that we are in an idealized setting.
Adaptive techniques, where the levels are computed on-the-fly have been introduced in the
2000s in various contexts, under different names: Adaptive Multilevel Splitting (AMS) [5], [6],
[7], Subset simulation [2] and Nested sampling [13].
In this paper, we focus on the versions of AMS algorithms studied in [3]. Two parameters
are required: a number of (interacting) replicas n, and a fixed integer k∈{1,...,n−1}, such
that the proportion of replicas that are killed and resampled at each iteration is k/n. The
version with k=1 has been studied in [10], and is also (in the idealized setting) a special case
of the Adaptive Last Particle Algorithm of [14].
Afamilyofestimators(pˆn,k) isintroducedin[3]. Themainpropertyestablished
n≥2,1≤k≤n−1
there is unbiasedness: E[pˆn,k]=p. Moreover, the computational cost is analyzed in the regime
when k is fixed and n→+∞. Nevertheless comparisons between different values of k is made
assuming a non natural procedure: N independent realizations of the algorithm are necessary
todefineanempiricalestimatorwithrealizationsofpˆn,k. Weremovethisprocedurebyshowing
directly an asymptotic normality result for the estimator pˆn,k, allowing for the direct use of
Gaussian asymptotic confidence intervals.
Other Central Limit Theorems for Adaptive Multilevel Splitting estimators (in different
parameter regimes) have been obtained in [4], [5] and [8].
The main result of this paper is Theorem 2: if k and a are fixed, under the assumption
that the cumulative distribution function of X is continuous, when n → +∞ we have the
convergence in law of √n(cid:0)pˆn,k−p(cid:1) to a centered Gaussian random variable, with variance
−p2log(p) (independent on k).
Themainnoveltyofthepaperisthecasek>1: indeedwhenk=1thelawoftheestimator
is explicitly known (it involves a Poisson random variable with parameter −nlog(p)) and the
Central Limit Theorem can be derived by hand. More precisely, we prove the asymptotic
normality of log(pˆn,k), and conclude thanks to the delta-method. However when k >1, the
key idea is to prove a functional equation, as introduced in [3], for the characteristic function
of this random variable; the basic ingredient is a decomposition according to the first step of
the algorithm. Thus one of the main messages of this paper is that the functional equation
technique allows to prove several key properties of the AMS algorithm in the idealized setting:
unbiasedness and asymptotic normality.
Thepaperisorganizedasfollows. InSection2,weintroducethemainobjects: theidealized
setting (Section 2.1) and the AMS algorithm (Section 2.2). Our main result (Theorem 2) is
statedinSection2.3. Section3isdevotedtothedetailedproofofthisresult. FinallySection4
contains a numerical illustration of the Theorem.
2 Adaptive Multilevel Splitting Algorithms
2.1 Setting
Let X be some real random variable. We assume that X > 0 almost surely. We want to
estimate the probability p=P(X >a), where a>0 is some threshold. When a goes to +∞, p
goes to 0 and we have to estimate a rare event. More generally, we introduce the conditional
probability for 0≤x≤a
P(x)=P(X >a|X >x). (1)
We notice that p=P(0) and that P(a)=1.
The following assumption is crucial:
2
Assumption 1 Let F denote the cumulative distribution function of X. We assume that F
is continuous.
2.2 Algorithm
The algorithm depends on two parameters, fixed once and for all:
• the number of replicas n;
• the number k∈{1,...,n−1} of replicas that are resampled at each iteration.
The other necessary parameters are the initial condition x and the stopping threshold a.
In the sequel, when we consider a random variable Xj, the subscript i denotes the index in
i
{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 ,...,Y ) independent and identically distributed (i.i.d.) real valued
1 n
random variables with continuous cumulative distribution function, there exists almost surely
a unique (random) permutation σ of {1,...,n} such that Y < ... < Y . For any
σ(1) σ(n)
k ∈ {1,...,n}, we then use the classical notation Y = Y to denote the k-th order
(k) σ(k)
statistics of the sample Y.
Algorithm 1 (Adaptive Multilevel Splitting)
Initialization: DefineZ0 =x. Sampleni.i.d. realizationsX0,...,X0,withthelawL(X|X >
1 n
x).
Define Z1 =X0 , the k-th order statistics of the sample X0 =(X0,...,X0), and σ1 the
(k) 1 n
(a.s.) unique associated permutation: X0 <...<X0 .
σ1(1) σ1(n)
Set j =1.
Iterations (on j ≥1): While Zj <a:
• Conditionally on Zj, sample k new independent random variables (Yj,...,Yj), according
1 k
to the law L(X|X >Zj).
• Set
(cid:40)
Yj if (σj)−1(i)≤k
Xj = (σj)−1(i)
i Xj−1 if (σj)−1(i)>k.
i
In other words, the particle with index i is killed and resampled according to the law
L(X|X > Zj) if Xj−1 ≤ Zj, and remains unchanged if Xj−1 > Zj. Notice that the
condition (σj)−1(i)i≤k is equivalent to i∈(cid:8)σj(1),...,σj(k)i(cid:9).
• Define Zj+1 =Xj , the k-th order statistics of the sample Xj =(Xj,...,Xj), and σj+1
(k) 1 n
the (a.s.) unique associated permutation: Xj <...<Xj .
σj+1(1) σj+1(n)
• Finally increment j ←j+1.
End of the algorithm: Define Jn,k(x)=j−1 as the (random) number of iterations. Notice
that Jn,k(x) is such that ZJn,k(x) <a and ZJn,k(x)+1 ≥a.
The estimator of the probability P(x) is defined by
(cid:18) k(cid:19)Jn,k(x)
pˆn,k(x)=Cn,k(x) 1− , (2)
n
with
Cn,k(x)= 1Card(cid:110)i; XJn,k(x) ≥a(cid:111). (3)
n i
When x=0, to simplify notations we set pˆn,k =pˆn,k(0).
2.3 The Central Limit Theorem
The main result of the paper is the following asymptotic normality result.
Theorem 2 Under assumption 1, we have the following convergence in law, when n→+∞,
for any fixed k and a:
√
n(cid:0)pˆn,k−p(cid:1)→N(cid:0)0,−p2log(p)(cid:1). (4)
3
Noticethattheasymptoticvariancedoesnotdependonk. Thisresultcanbeusedtodefine
asymptotic normal confidence intervals, for one realization of the algorithm and n → +∞.
However, the speed of convergence is not known and may depend on the estimated probability
n and on the parameter k.
Thanks to Theorem 2, we can study the cost of the use of one realization of the AMS
algorithm to obtain a given accuracy when n→+∞. In [3], the cost was analyzed when using
a sample of M independent realizations of the algorithm, giving an empirical estimator, and
the analysis was based on an asymptotic analysis of the variance in the large n limit.
Let(cid:15)besomefixedtoleranceerror,andα>0. Denoter suchthatP(Z ∈[−r ,r ])=1−α,
α α α
where Z is a standard Gaussian random variable.
Then for n large, an asymptotic confidence interval with level 1−α, centered around p, is
(cid:34) (cid:112) (cid:112) (cid:35)
r −p2log(p) r −p2log(p)
p− α √ ,p+ α √ .
n n
Then the (cid:15) error criterion |pˆn,k−p|≤(cid:15) is achieved for n of size −p2log(p)rα2.
(cid:15)2
However,inaverageonerealizationoftheAMSalgorithmrequiresanumberofstepsofthe
order −nlog(p)/k, with k random variables sampled at each iteration. Another cost comes
fromthesortingofthereplicasatinitialization,andtheinsertionateachiterationofthek new
sampled replicas in the sorted ensemble of the non-resampled ones. Thus the cost to achieve
an accuracy of size (cid:15) is in the large n regime of size nlog(n)(cid:0)−p2log(p)(cid:1), not depending on k.
This cost must be compared with the one when using a pure Monte-Carlo approximation,
with an ensemble of non-interacting replicas of size n: thanks to the Central Limit Theorem,
the tolerance criterion error (cid:15) is satisfied for n of size −p(1−p)rα2. Despite the log(n) factor in
(cid:15)2
the AMS case, the performance is improved since p2log(p)=o(p) when p→0.
Remark 1 Unlike in [3], here we are not able to compare the different choices for k, since at
the limit n→+∞ all the involved quantities do not depend on k.
3 Proof of the Central Limit Theorem
The proof is divided into the following steps. First, thanks to Assumption 1, we reduce our
study to the case when X is distributed according to the exponential law with parameter
1: P(X > z) = exp(−z) for any z > 0. Second, we explain why it is simpler to look at
log(pˆn,k(x))insteadofpˆn,k(x),andtousethedelta-methodtogettheresult. Thethirdstepis
the introduction of the characteristic function of log(pˆn,k(x)); following the definition of the
algorithm, we prove that it is solution of a functional equation with respect to x, which can be
transformed into a linear ODE of order k. Finally, we analyze the solution of this ODE in the
limit n→+∞.
3.1 Reduction to the exponential case
WestartwithareductionofthestudytothecasewhentherandomvariableX isexponentially
distributed with parameter 1. It is based on a change of variable with the following function:
(cid:0) (cid:1)
Λ(x)=−log 1−F(x) . (5)
It is well-known that F(X) is distributed according to the uniform law on (0,1) (thanks to the
continuity Assumption 1), and thus Λ(X) is exponentially distributed with parameter 1. In
Section 3 (see Corollary 3.4) of [3], it is proved that, the law of the estimator pˆn,k is equal to
the one of qˆn,k, which is the estimator defined with the same values of the parameters n and
k, but with two differences: the random variable is exponentially distributed with parameter
1, and the stopping level is Λ(a). Notice that E(cid:2)qˆn,k(cid:3)=exp(cid:0)−Λ(a)(cid:1)=1−F(a)=p (by the
unbiasedness result of [3]).
Since the arguments are intricate, we do not repeat them here and we refer the interested
reader to [3]; from now on, we thus assume:
Assumption 3 Assume that X is exponentially distributed with parameter 1: we denote
L(X)=E(1).
As a consequence, it is enough to prove the following result:
4
Proposition 1 When n→+∞
√
n(cid:0)pˆn,k−p(cid:1)→N(cid:0)0,aexp(−2a)(cid:1). (6)
We will need the following notations:
(cid:0) (cid:1)
• f(z) = exp(−z)1 (resp. F(z) = 1−exp(−z) 1 ) for the density (resp. the
z>0 z>0
cumulative distribution function) of the exponential law E(1) with parameter 1.
• f (z)=k(cid:0)n(cid:1)F(z)k−1f(z)(cid:0)1−F(z)(cid:1)n−k for the density of the k-th order statistics of a
n,k k
sample of size n, made of independent and exponentially distributed random variables,
with parameter 1.
Finally, in order to deal with the conditional distributions L(X|X >x) (under Assumption
3itisjustashiftedexponentialdistributionx+E(1)thankstothelossofmemorypropertyof
the exponential law) in the algorithm, we set for any x≥0 and any y≥0
f(y;x)=f(y−x), F(y;x)=F(y−x),
f (y;x)=f (y−x),
n,k n,k (7)
(cid:90) y
F (y)= f (z)dz, F (y;x)=F (y−x).
n,k n,k n,k n,k
−∞
The following result will be very useful:
 d
 dxfn,1(y;x)=nfn,1(y;x).
(8)
d
 for k∈{2,...,n−1}, dxfn,k(y;x)=(n−k+1)(fn,k(y;x)−fn,k−1(y;x)).
The proof follows from straightforward computations.
3.2 Proof of the CLT
The first important idea is to prove the Proposition 1 for any initial condition x∈[0,a], not
only for x=0:
√
n(cid:0)pˆn,k(x)−p(x)(cid:1)→N(cid:0)0,(a−x)exp(−2(a−x))(cid:1). (9)
A natural idea is to introduce the characteristic function of pˆn,k(x), and to follow the
strategy developed in the previous work. Nevertheless, we are not able to derive a useful
functionalequationwithrespecttothexvariable. Abetterstrategyistolookatthelogarithm
of the estimator, and to use a particular case of the delta-method (see for instance [15],
S√enc(cid:0)tiθon−3θ):)→if fNor(0so,mσ2e),stehqeunen√cne(cid:0)oefxrpe(aθl)ra−ndexopm(θv)a(cid:1)r→iabNle(cid:0)s0(,θenx)pn(∈2Nθ)aσn2d(cid:1)-sobmotehθin∈dRistwriebuhtaiovne
n n
when n→+∞.
We thus introduce for any t∈R and any 0≤x≤a
φ (t,x):=E(cid:104)exp(cid:16)it√n(cid:0)log(pˆn,k(x))−log(p(x))(cid:1)(cid:17)(cid:105). (10)
n,k
For technical reasons, we also set
χ (t,x):=E(cid:104)exp(cid:16)it√npˆn,k(x)(cid:17)(cid:105)=exp(cid:0)it√n(x−a)(cid:1)φ (t,x), (11)
n,k n,k
on which we derive a functional equation, in order to prove the following result, which implies
Proposition 1 with the choice x=0.
Proposition 2 When n→+∞, for any 0≤x≤a and any t∈R
(cid:18)t2(x−a)(cid:19)
φ (t,x)→exp . (12)
n,k 2
The rest of this section is devoted to the proof of Proposition 2, thanks to the following
series of results.
5
Lemma 1 (Functional Equation) For any n∈N and any k∈{1,...,n−1}, and for any
t∈R,thefunctionx(cid:55)→χ (t,x)issolutionofthefollowingfunctionalequation(withunknown
n,k
χ): for any 0≤x≤a
χ(t,x) = eit√nlog(1−nk)(cid:90) aχ(t,y)fn,k(y;x)dy (13)
x
+k(cid:88)−1eit√nlog(1−nl)P(S(x)n(l) <a≤S(x)n(l+1)), (14)
l=0
where (S(x)n) are iid with law L(X|X >x) and where S(x)n is the l-th order statistics
j 1≤j≤n (l)
of this sample (with convention S(x)n =x).
(0)
Remark 2 The function (cid:55)→φ (t,x) is solution of a more complex functional equation: for
n,k
any 0≤x≤a
(cid:90) a √ √
φ (t;x) = eit nlog(1−k/n)e−it n(x−y)φ (t;y)f (y;x)dy
n,k n,k n,k
x
+k(cid:88)−1e−it√n(x−a)eit√nlog(1−nl)P(S(x)n(l) <a≤S(x)n(l+1)).
l=0
The advantage of the equation on χ is that we are able to deduce that x (cid:55)→ χ (t,x) is
n,k n,k
solution of a linear ODE with constant coefficients. On x(cid:55)→φ (t,x) the ODE would not have
n,k
constant coefficients and would be more difficult to express.
Remark 3 If instead of considering log(pˆn,k)(x) we try to prove directly the central limit
theorem on pˆn,k(x), the same strategy leads to a functional equation where on the left-hand side
we have the parameter t, and on the right-hand side, in the integrand we have the parameter
(cid:0) (cid:1)
t 1−k/n .
Proof. The idea (like in the proof of Proposition 4.2 in [3]) is to decompose the expectation
according to the value of the first level Z1 =X0 . On the event (cid:8)Z1 >a(cid:9)=(cid:8)Jn,k(x)=0(cid:9),
(k)
the algorithm stops and pˆn,k(x)= n−l for the unique l ∈{0,...,k−1} such that S(x)n <
n (l)
a≤S(x)n . Thus
(l+1)
E[eit√nlog(pˆn,k(x))1Jn,k(x)=0]=k(cid:88)−1eit√nlog(1−nl)P(S(x)n(l) <a≤S(x)n(l+1)). (15)
l=0
If Z1 <a, for the next iteration the algorithm restarts at the initial condition Z1, and
E[eit√nlog(pˆn,k(x))1 ]
Jn,k(x)>0
=E(cid:34)eit√nlog(1−nk)E[eit√nlog(cid:18)Cn,k(x)(1−nk)Jn,k(x)−1(cid:19)|Z1]1Z1<a(cid:35)
=eit√nlog(1−nk)E(cid:104)E[eit√nlog(pˆn,k(Z1))|Z1]1Z1<a(cid:105) (16)
=eit√nlog(1−nk)E(cid:2)χn,k(t,Z1)1Z1<a(cid:3)
=eit√nlog(1−nk)(cid:90) aχn,k(t,y)fn,k(y;x)dy.
x
With the two identities and the definition (11) of χ , we thus obtain (13).
n,k
(cid:3)
We exploit the functional equation (13) on x(cid:55)→χ (t,x), to prove that this function is
n,k
solution of a linear ODE with constant coefficients.
Lemma 2 (ODE) Let n and k ∈ {1,...,n−2} be fixed. There exist real numbers µn,k
and (rn,k) , depending only on n and k, such that for any fixed t ∈ R, the function
m 0≤m≤k−1
x(cid:55)→χ (t,x) satisfy the following Linear Ordinary Differential Equation (ODE) of order k:
n,k
for x∈[0,a]
ddxkkχn,k(t,x)=eit√nlog(1−nk)µn,kχn,k(t,x)+ k(cid:88)−1rmn,kddxmmχn,k(t,x). (17)
m=0
6
The coefficients µn,k and (rn,k) satisfy the following properties:
m 0≤m≤k−1
µn,k =(−1)kn...(n−k+1)
k−1 (18)
λk− (cid:88)rn,kλm =(λ−n)...(λ−n+k−1) for all λ∈R.
m
m=0
Proof. The proof follows the same lines as Proposition 6.4 in [3]. We introduce the following
notation
Θn,k(t,x):=k(cid:88)−1eit√nlog(1−nl)P(S(x)n(l) <a≤S(x)n(l+1)).
l=0
Then by recursion, we can prove that for 0≤l≤k−1 and for any x≤a and t∈R
ddxll (χn,k(t,x)−Θn,k(t,x)) = µnl,keit√nlog(1−nk)(cid:90) aχn,k(t,y)fn,k−l(y;x)dy
x
+(cid:88)l−1 rn,k dm (χ (t,x)−Θ (t,x)). (19)
m,ldxm n,k n,k
m=0
The idea for the recursion is to differentiate (19) with respect to x and to use (8). The
coefficients satisfy the following recursion (notice that they do not depend on t):
µn,k =1,µn,k =−(n−k+l+1)µn,k;
0 l+1 l
rn,k =−(n−k+l+1)rn,k, if l>0,
 0,l+1 0,l (20)
rn,k =rn,k −(n−k+l+1)rn,k, 1≤m≤l,
m,l+1 m−1,l m,l
rn,k =−1.
l,l
Differentiating once more when l=k−1, we obtain
ddxkk (χn,k(t,x)−Θn,k(t,x)) = µn,keit√nlog(1−nk)χn,k(t,x)
+k(cid:88)−1rn,k dm (χ (t,x)−Θ (t,x)), (21)
m dxm n,k n,k
m=0
with µn,k :=µn,k and rn,k :=rn,k. These coefficients are the same as in [3]; in particular it
k m m,k
was proved there that the claimed properties are satisfied (using the unbiasedness result for
the estimation of the probability, see Section 6.4 there). The polynomial equality of (18) is in
fact equivalent to the following identity: for all j ∈{0,...,k−1}
dk exp((n−k+j+1)(x−a))= k(cid:88)−1rn,k dm exp((n−k+j+1)(x−a)).
dxk m dxm
m=0
Using the expression of the order statistics in the exponential case, Θ (t,.) is a linear
n,k
combination of the exponential functions z(cid:55)→exp(−nz),...,exp(−(n−k+1)z); therefore
dk Θ (t,x)= k(cid:88)−1rn,k dm Θ (t,x),
dxk n,k m dxm n,k
m=0
and thus (21) gives (17).
(cid:3)
To conclude, it remains to express the solution of (17), and to analyze the asymptotics
n → +∞. Since the ODE is of order k, in order to uniquely determine the solution, more
information is required: for instance we need to know the vector of the derivatives of order
0,1,...,k−1 of x(cid:55)→χ (t,x) at some point. From the functional equation (13), we are able
n,k
to show the following asymptotic result at point a, when n→+∞.
Lemma 3 (Initial condition) For any fixed k∈{1,...,} and any t∈R, we have
(cid:40)
χ (t,a)=1
n,k
ddxmmχn,k(t,x)(cid:12)(cid:12)(cid:12)x=a n→=∞O(√1n)nm if m∈{1,...,k−1}. (22)
7
Proof. The equality χ (t,a) = 1 is trivial, since pˆn,k(a) = 1. From (19) and (21), we
n,k
immediately get (by recursion) that for 1≤m≤k−1
dm (cid:12) dm (cid:12)
χ (t,x)(cid:12) = Θ (t,x)(cid:12) .
dxm n,k (cid:12)x=a dxm n,k (cid:12)x=a
We introduce the following decomposition
Θn,k(t,x)=k(cid:88)−1(eit√nlog(1−nl)P(S(x)n(l) <a≤S(x)n(l+1))
l=0
=k(cid:88)−1(cid:16)eit√nlog(1−nl)−1(cid:17)(cid:0)Fn,l(a;x)−Fn,l+1(a;x)(cid:1)
l=0
k−1
+(cid:88)P(S(x)n <a≤S(x)n )
(l) (l+1)
l=0
=:Ω (t,x)+1−F (a;x),
n,k n,k
where F denotes the cumulative distribution function of the l-th order statistics (with the
n,l
convention F (a;x)=1 for x≤a).
n,0
Thanks to (8) and a simple recursion on l, it easy to prove that for any 0≤l≤k and any
m≥1
dm (cid:12)
F (a;x)(cid:12) =O(nm), (23)
dxm n,l (cid:12)x=a
which immediately yields
dm (cid:12) 1
Ω (t,x)(cid:12) = O(√ )nm.
dxm n,k (cid:12)x=a n→∞ n
In fact, it is possible to prove a stronger result: if 1≤l≤k and 0≤m<l then
dm (cid:12)
F (a;x)(cid:12) =0,
dxm n,l (cid:12)x=a
by recursion on l and using (8) recursively on m. We thus obtain for 1≤m≤k−1
dm (cid:0)1−F (a;x)(cid:1)(cid:12)(cid:12) =0.
dxm n,k (cid:12)x=a
This concludes the proof of the Lemma.
(cid:3)
The last useful result is the following:
Lemma 4 (Asymptotic expansion) Let k∈{1,...,} and t∈R be fixed. Then for n large
enough, we have
k
χn,k(t,x)=(cid:88)ηnl,k(t)eλln,k(t)(x−a), (24)
l=1
for complex coefficients satisfying:
 √
 λ1 (t) = it n+ t2 +o(1),
n,k n→∞ 2 (25)
η1 (t) → 1;
 n,k
n→∞
and for 2≤l≤k

 λln,k(t) ∼ n(1−ei2π(kl−1)),
n→∞ (26)
 ηnl,k(t) → 0.
n→∞
Proof. We denote by (λl (t)) the roots of the characteristic equation (with unknown
n,k 1≤l≤k
λ∈C):
(n−λ)...(n−k+1−λ) −eit√nlog(1−nk) =0
n...(n−k+1)
Bythecontinuitypropertyoftherootsofacomplexpolynomialofdegreek withrespecttoits
coefficients, we have
λl (t)
λl (t):= n,k → λl ,
n,k n n→∞ ∞
8
where (λl (t)) are the roots of (1−λ)k =1: thus, we get λ1 (t) = o(n),
∞ 1≤l≤k n,k
n→∞
λln,k(t) ∼ n(1−ei2π(kl−1)).
n→∞
To study more precisely the asymptotic behavior of λ1 (t), we postulate an ansatz
n,k
√
λ1 (t) = c n+d +o(1),
n,k t t
n→∞
and identify the coefficients c =it and d =t2/2 thanks to
t t
(cid:18)1− √ct − dt +o(1)(cid:19)k = 1− √ctk − dtk−(cid:0)k2(cid:1)c2t +o(cid:18)1(cid:19)
n n n→∞ n n n
eit√nlog(1−nk) = 1− √itk − t2k2 +o(cid:18)1(cid:19).
n→∞ n 2n n
In particular, for n large enough, (λl (t)) are pairwise distinct, and (24) follows.
n,k 1≤l≤k
The coefficients (ηl (t)) are then solutions of the following linear system of order k:
n,k 1≤l≤k
 η1 (t)+...+ηk (t)=χ (a),
 ηnn1,,kk(t)λ1n,k(t)+n.,.k.+ηnk,k(nt,)kλkn,k(t)= n1ddxχn,k(t,a), (27)
...
 η1 (t)(cid:16)λ1 (t)(cid:17)k−1+...+ηk (t)(cid:16)λk (t)(cid:17)k−1 = 1 dk−1 χ (t,a).
n,k n,k n,k n,k nk−1dxk−1 n,k
Using Cramer’s rule, we express each ηl (t) as a ratio of determinants (the denominator
n,k
being a Vandermonde determinant is non zero when n is large enough). For l∈{2,...,k}, we
have
det(Ml (t))
ηl (t)= n,k → 0
n,k V(cid:16)λ1 (t),...,λk (t)(cid:17) n→+∞
n,k n,k
where
 1 1 ... 1 ... 1 
Mnl,k(t)= λ1n,k...(t) λ2n,k...(t) ...... O(√...1n) ...... λkn,k...(t) 
 
(cid:0)λ1n,k(t)(cid:1)k−1 (cid:0)λ2n,k(t)(cid:1)k−1 ... O(√1n) ... (cid:0)λkn,k(t)(cid:1)k−1
is such that det(Ml (t)) → 0 (since λ1 (t) → 0), and the denominator is given by a
n,k n,k
n→+∞
(cid:16) (cid:17) (cid:16) (cid:17)
1 k 1 k
Vandermonde determinant V λ (t),...,λ (t) → V λ (t),...,λ (t) (cid:54)=0.
n,k n,k ∞ ∞
n→+∞
Finally, η1 (t)=1−(cid:80)k ηl (t) → 1.
n,k l=2 n,k n→+∞
(cid:3)
The proof of Proposition 2 is now straightforward, recalling from (10) and (11) the relation
(cid:0) √ (cid:1)
φ (t,x) = exp −it n(x−a) χ (t,x) between the characteristic function φ and the
n,k n,k n,k
auxiliary function χ , and taking the limit n→+∞ thanks to Lemma 4.
n,k
4 Numerical results
In this Section, we provide some numerical illustrations of the Central Limit Theorem 2. We
apply the algorithm in the idealized setting with an exponentially distributed random variable
with parameter 1; indeed, we have seen that in the idealized setting, in law the algorithm
behaves like this, if the threshold is suitably chosen.
Inthesimulationswepresent,theprobabilityweestimateisexp(−6),whichisapproximately
2.10−3.
In figure 1, we fix the value k = 10, and we show histograms for n = 102,103,104,
with different values for the number M independent realizations of the algorithm, such that
nM =108 (we thus have empirical variance of the same order for all cases). When n increases,
thenormalityoftheestimatorisconfirmed. Infigure1,wegiveQ-Qplots,wheretheempirical
9
Figure 1: Histograms for k =10 and p=exp(−6): n=102,103,104 from left to right
Figure 2: Q-Q plot for k =10 and p=exp(−6): n=102,103,104 from left to right
quantilesofthesamplearecomparedwiththeexactquantilesofthestandardGaussianrandom
variable (after normalization).
InFigure3,weshowhistogramsforM =104independentrealizationsoftheAMSalgorithm
with n = 104 and k ∈ {1,10,100}; we also provide associated Q-Q plots in figure 4, which
prove that in this regime the normality assumption seems to be satisfied for all values of k.
Figure 3: Histograms for n=104 and p=exp(−6): k =1,10,100 from left to right
Figure 4: Q-Q plot for n=104 and p=exp(−6): k =1,10,100 from left to right
Acknowledgements: C.-E.B.wouldliketothankG.Samaey,T.Leli`evreandM.Rousset
for theinvitationto givea talkon thetopic ofthis paperat the11-th MCQMCConference, in
the special session on Mathematical aspects of Monte Carlo methods for molecular dynamics.
10