2.1. Loi des grands nombres, Théorème Central Limit#

On rappelle les deux théorèmes principaux qui justifient la méthode de Monte Carlo. Soit (Ω,A,P) un espace de probabilité.

Theorem 2.1 (Loi Forte des Grands Nombres, LFGN)

Soit (Xn)n1 une suite de variable aléatoire i.i.d. intégrable, i.e. E[|X1|]<+ alors

(2.1)#limn1nk=1nXk=E[X1]p.s.et dans L1.

La méthode de Monte Carlo est une application directe de cette loi des grands nombres. On appelle estimateur de Monte Carlo un estimateur In de la forme In=1nk=1nXk. Voyons une application en intégration numérique.

Soit f une fonction mesurable telle que 01|f(x)|dx<+. Pour obtenir une approximation de I=01f(x)dx, on écrit cette quantité sous forme probabiliste I=E[f(U)] avec UU([0,1]). On considère alors une suite (Un)n1 i.i.d. uniforme sur [0,1] et l’estimateur Monte Carlo

In=1nk=1nf(Uk)vérifielimnIn=Ip.s.

Plus généralement si on veut calculer I=Rdf(x)g(x)dxg est une densité de probabilité g:RdR+ et f:RdR telle que Rd|f(x)|g(x)dx<+, alors on écrit I=E[f(Y)] avec Y variable aléatoire de densité g (on note Yg), et on considère l’estimateur Monte Carlo

In=1nk=1nf(Yk)avec (Yn)n1 suite i.i.d. de loi g.

En pratique il faut arrêter la somme à une itération n donnée. Il faut donc connaitre l’erreur 1nk=1nXkE[X1]. Le deuxième théorème fondamental à connaitre est celui qui donne le comportement de cette erreur lorsqu’on la renormalise par la bonne quantité n.

Theorem 2.2 (Théorème Central Limite, TCL)

Soit (Xn)n1 une suite de variables aléatoires réelles i.i.d. de carré intégrable i.e. E[X2]<+. On note m=E[X1] et σ2=var(X1). Alors

(2.2)#n(1nk=1nXkm)LN(0,σ2).

Deux informations importantes apparaissent avec ce théorème:

  • Vitesse. La vitesse de convergence de In vers m est en n ce qui signifie qu’il faut multiplier par 100 le nombre d’itérations pour diviser l’erreur par 10. C’est une vitesse très lente comparer à d’autres méthodes numériques mais le principal avantage est que cette vitesse ne dépend pas de la dimension du problème (à la différence des méthodes classiques d’intégration numérique).

  • Variance. Plus la variance est grande, moins la méthode est efficace. Il faut donc construire des estimateurs Monte Carlo avec la variance la plus petite possible. Le rôle de la variance est primordiale et apparait dans le “critère d’arrêt” (le choix de n) basé sur les intervalles des confiance. Il est obligatoire de calculer la variance d’une méthode de Monte Carlo.

2.2. Intervalles de confiance#

Voici la définition classique d’un intervalle de confiance que l’on peut trouver un ouvrage de statistiques.

Definition 2.1 (Intervalle de confiance)

Soit X1,,Xn un échantillon de loi Pθ et α]0,1[ (α petit). Un intervalle de confiance pour θ de niveau de confiance 1α est un intervalle qui dépend de l’échantillon tel que la probabilité que θ soit dans cet intervalle soit au moins égale à 1α.

On va déduire tu TCL la construction d’un intervalle de confiance pour la moyenne m=E[X1].

Dans un premier temps, supposons la variance σ2 connue. D’après le TCL on a a>0,

(2.3)#limnP(nσ(1nk=1nXkm)[a,a])=P(G[a,a])

avec GN(0,1). On en déduit un intervalle pour m qui dépend de X1,,Xn

limnP(m[1nk=1nXk±aσn])=P(G[a,a])

Pour avoir un niveau de confiance de 1α il faut prendre a>0 tel que P(G[a,a])=1α. On utilise la fonction de répartition FN(0,1) et son inverse FN(0,1)1 de sorte que

P(G[a,a])=1P(|G|>a)=12(1FN(0,1)(a))=1α

ce qui indique qu’il faut considérer

a=FN(0,1)1(1α2).

Dans la suite on note qα=FN(0,1)1(1α2) de sorte que

(2.4)#limnP(m[1nk=1nXk±qασn])=1α.

Pour aller plus loin il faut faire l’approximation suivante (justifiée par le théorème de Berry-Essen, cf. plus loin) pour n suffisamment grand,

(2.5)#P(m[1nk=1nXk±qασn])1α,

ce qui justifie l’intervalle de confiance (asymptotique) suivant

m[1nk=1nXk±qασn]avec probabilité 1α.

Lorsqu’on prend α=0.05 c’est à dire que l’on veut un niveau de confiance de 95% on a qα=1.96.

Si la variance σ2 est inconne, on utilise un estimateur Monte Carlo pour estimer σ2=var(X1). Par exemple on peut considérer

(2.6)#σn2=1n1k=1n(XkIn)2avecIn=1nk=1nXk,=1n1k=1nXk2nn1In2,

(il s’agit de l’estimateur sans biais de la variance). Ainsi par la LFGN limnσn2=σ2 p.s., et par le lemme de Slutsky on en déduit le TCL suivant

nσn(1nk=1nXkm)LN(0,1).

Comme précédemment on peut en déduire l’intervalle de confiance (asymptotique) suivant

m[1nk=1nXk±qασnn]avec probabilité 1α.

Revenons à l’approximation (2.5) qui est plus forte que juste l’application du TCL. En fait cette approximation est justifée par un résultat non-asymptotique (vrai pour tout n1) appelé théorème de Berry-Essen.

Theorem 2.3 (Berry-Essen)

Soit (Xn)n1 une suite de variables aléatoires rélles i.i.d. telle que E[X1]=0, σ2=var(X1) et ρ=E[|X1|3]<+. Soit Fn la fonction de répartition de l’erreur renormalisée 1σnk=1nXk. Alors

xR,n1,|Fn(x)FN(0,1)(x)|3ρσ3n.

2.2.1. IC non asymptotique via Tchebychev#

Pour compléter cette partie sur les intervalles de confiances (I.C.), comparons l’I.C. asymptotique obtenu par le TCL avec d’autres I.C. Prenons par exemple l’inégalité de Tchebychev qui permet de construire un I.C. non asymptotique.

Proposition 2.1 (Tchebychev)

Si X est de carré intégrable, alors

r>0,P(|XE[X]|r)var(X)r2.

Appliquons cette proposition avec X=1nk=1nXk. Alors pour tout r>0

P(|1nk=1nXkm|r)var(1nk=1nXk)r2=σ2nr2,

et

n1,r>0,P(m[1nk=1nXk±r])1σ2nr2.

Pour en déduire un I.C. de niveau 1α il faut prendre r tel que σ2nr2=α c’est à dire r=1/ασn. Si on choisit α=0.05 la constante 1/α qui apparait devant σn est égale à 4.47. Pour le TCL on rappelle que la constante devant σn vaut qα=1.96.

2.2.2. IC non asymptotique via Chernov#

Si on suppose que X a des moments exponentiels, i.e. τ>0, E[eτX]<+ il est possible d’obtenir un I.C. non asymptotique plus reserré en considérant par exemple l’inégalité de Chernov.

Proposition 2.2 (Chernov)

Soit (Xn)n1 une suite de variables aléatoires réelles i.i.d. de loi μ ayant des moments exponentiels et m=E[X1]. Alors

a>m,P(1nk=1nXk>a)enhμ(a),a<m,P(1nk=1nXk<a)enhμ(a),

hμ est la transformée de Cramer de la loi μ définie par hμ(x)=supθR(θxψ(θ)) et ψ(θ)=logE[eθX].

Appliquons cette proposition dans le cas Gaussien uniquement. On suppose donc que X1N(0,σ2) (on prend m=0 pour simplifier). D’après cette proposition on a pour tout ε>0,

(2.7)#P(|1nk=1nXk|>ε)=P(1nk=1nXk>ε)+P(1nk=1nXk<ε)2enhμ(ε).

La transformée de Cramer de la loi N(0,σ2) est donnée par hμ(x)=12(xσ)2 donc

P(|1nk=1nXk|ε)12en2(εσ)2.

Pour construire un I.C. non asymptotique de niveau 1α on choisit alors ε=2log(α/2)σn. La constante pour α=0.05 vaut alors 2log(0.025)=2.716 (à comparer avec 1.96 et 4.47).