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 \((\Omega, \mathcal{A}, \mathbf{P})\) un espace de probabilité.

Theorem 2.1 (Loi Forte des Grands Nombres, LFGN)

Soit \((X_n)_{n \ge 1}\) une suite de variable aléatoire i.i.d. intégrable, i.e. \(\mathbf{E}[|X_1|] < +\infty\) alors

(2.1)\[\begin{equation} \lim_{n \to \infty} \frac{1}{n} \sum_{k = 1}^n X_k = \mathbf{E}[X_1] \quad p.s. \quad \text{et dans $\mathbf{L}^1$.} \end{equation}\]

La méthode de Monte Carlo est une application directe de cette loi des grands nombres. On appelle estimateur de Monte Carlo un estimateur \(I_n\) de la forme \(I_n = \frac{1}{n} \sum_{k=1}^n X_k\). Voyons une application en intégration numérique.

Soit \(f\) une fonction mesurable telle que \(\int_0^1 |f(x)| \operatorname{d}\!x < +\infty\). Pour obtenir une approximation de \(I = \int_0^1 f(x) \operatorname{d}\!x\), on écrit cette quantité sous forme probabiliste \(I = \mathbf{E}[f(U)]\) avec \(U \sim \mathcal{U}([0,1])\). On considère alors une suite \((U_n)_{n \ge 1}\) i.i.d. uniforme sur \([0,1]\) et l’estimateur Monte Carlo

\[\begin{equation*} I_n = \frac{1}{n} \sum_{k=1}^n f(U_k) \quad \text{vérifie} \quad \lim_n I_n = I \quad p.s. \end{equation*}\]

Plus généralement si on veut calculer \(I = \int_{\mathbf{R}^d} f(x) g(x) \operatorname{d}\!x\) où \(g\) est une densité de probabilité \(g:\mathbf{R}^d \to \mathbf{R}_+\) et \(f:\mathbf{R}^d \to \mathbf{R}\) telle que \(\int_{\mathbf{R}^d} |f(x)| g(x) \operatorname{d}\!x < +\infty\), alors on écrit \(I = \mathbf{E}[f(Y)]\) avec \(Y\) variable aléatoire de densité \(g\) (on note \(Y \sim g\)), et on considère l’estimateur Monte Carlo

\[\begin{equation*} I_n = \frac{1}{n} \sum_{k=1}^n f(Y_k) \quad \text{avec $(Y_n)_{n \ge 1}$ suite i.i.d. de loi $g$}. \end{equation*}\]

En pratique il faut arrêter la somme à une itération \(n\) donnée. Il faut donc connaitre l’erreur \(\frac{1}{n} \sum_{k=1}^n X_k - \mathbf{E}[X_1]\). 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é \(\sqrt{n}\).

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

Soit \((X_n)_{n \ge 1}\) une suite de variables aléatoires réelles i.i.d. de carré intégrable i.e. \(\mathbf{E}[X^2] < +\infty\). On note \(m = \mathbf{E}[X_1]\) et \(\sigma^2 = \operatorname{var}(X_1)\). Alors

(2.2)\[\begin{equation} \sqrt{n} \Big(\frac{1}{n} \sum_{k=1}^n X_k - m \Big) \xrightarrow{\mathcal{L}} \mathcal{N}(0, \sigma^2). \end{equation}\]

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

• Vitesse. La vitesse de convergence de \(I_n\) vers \(m\) est en \(\sqrt{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 \(X_1, \dots, X_n\) un échantillon de loi \(\mathbf{P}_\theta\) et \(\alpha \in ]0,1[\) (\(\alpha\) petit). Un intervalle de confiance pour \(\theta\) de niveau de confiance \(1-\alpha\) est un intervalle qui dépend de l’échantillon tel que la probabilité que \(\theta\) soit dans cet intervalle soit au moins égale à \(1-\alpha\).

On va déduire tu TCL la construction d’un intervalle de confiance pour la moyenne \(m = \mathbf{E}[X_1]\).

Dans un premier temps, supposons la variance \(\sigma^2\) connue. D’après le TCL on a \(\forall a > 0\),

(2.3)\[\begin{equation} \lim_{n \to \infty} \mathbf{P} \Bigl(\frac{\sqrt{n}}{\sigma} \big( \frac{1}{n} \sum_{k=1}^n X_k - m \big) \in [-a, a] \Bigr) = \mathbf{P}\big(G \in [-a, a]\big) \end{equation}\]

avec \(G \sim \mathcal{N}(0,1)\). On en déduit un intervalle pour \(m\) qui dépend de \(X_1, \dots, X_n\)

\[\begin{equation*} \lim_{n \to \infty} \mathbf{P} \Big(m \in \Big[\frac{1}{n} \sum_{k=1}^n X_k \pm \frac{a \sigma}{\sqrt{n}} \Big] \Big) = \mathbf{P}\big(G \in [-a, a]\big) \end{equation*}\]

Pour avoir un niveau de confiance de \(1-\alpha\) il faut prendre \(a > 0\) tel que \(\mathbf{P}\big( G \in [-a, a] \big) = 1 - \alpha\). On utilise la fonction de répartition \(F_{\mathcal{N}(0,1)}\) et son inverse \(F^{-1}_{\mathcal{N}(0,1)}\) de sorte que

\[\begin{equation*} \mathbf{P}\big( G \in [-a, a] \big) = 1 - \mathbf{P}\big(|G| > a\big) = 1 - 2 (1 - F_{\mathcal{N}(0,1)}(a)) = 1 - \alpha \end{equation*}\]

ce qui indique qu’il faut considérer

\[\begin{equation*} a = F_{\mathcal{N}(0,1)}^{-1}\big(1 - \frac{\alpha}{2} \big). \end{equation*}\]

Dans la suite on note \(q_\alpha = F_{\mathcal{N}(0,1)}^{-1}\big(1 - \frac{\alpha}{2} \big)\) de sorte que

(2.4)\[\begin{equation} \lim_{n \to \infty} \mathbf{P} \Big( m \in \Big[ \frac{1}{n} \sum_{k=1}^n X_k \pm \frac{q_\alpha \sigma}{\sqrt{n}} \Big] \Big) = 1 - \alpha. \end{equation}\]

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)\[ \mathbf{P} \Big( m \in \Big[ \frac{1}{n} \sum_{k=1}^n X_k \pm \frac{q_\alpha \sigma}{\sqrt{n}} \Big] \Big) \simeq 1 - \alpha,\]

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

\[\begin{equation*} m \in \Big[ \frac{1}{n} \sum_{k=1}^n X_k \pm \frac{q_\alpha \sigma}{\sqrt{n}} \Big] \quad \text{avec probabilité $1-\alpha$}. \end{equation*}\]

Lorsqu’on prend \(\alpha = 0.05\) c’est à dire que l’on veut un niveau de confiance de \(95\%\) on a \(q_\alpha = 1.96\).

Si la variance \(\sigma^2\) est inconne, on utilise un estimateur Monte Carlo pour estimer \(\sigma^2 = \operatorname{var}(X_1)\). Par exemple on peut considérer

(2.6)\[\begin{equation} \begin{split} \sigma_n^2 &= \frac{1}{n-1} \sum_{k=1}^n (X_k - I_n)^2 \quad \text{avec} \quad I_n = \frac{1}{n} \sum_{k=1}^n X_k, \\ & = \frac{1}{n-1} \sum_{k=1}^n X_k^2 - \frac{n}{n-1} I_n^2, \end{split} \end{equation}\]

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

\[\begin{equation*} \frac{\sqrt{n}}{\sigma_n} \Big( \frac{1}{n} \sum_{k=1}^n X_k - m \Big) \xrightarrow{\mathcal{L}} \mathcal{N}(0,1). \end{equation*}\]

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

\[\begin{equation*} m \in \Big[ \frac{1}{n} \sum_{k=1}^n X_k \pm \frac{q_\alpha \sigma_n}{\sqrt{n}} \Big] \quad \text{avec probabilité $1-\alpha$}. \end{equation*}\]

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 \(n \ge 1\)) appelé théorème de Berry-Essen.

Theorem 2.3 (Berry-Essen)

Soit \((X_n)_{n \ge 1}\) une suite de variables aléatoires rélles i.i.d. telle que \(\mathbf{E}[X_1] = 0\), \(\sigma^2 = \operatorname{var}(X_1)\) et \(\rho = \mathbf{E}[|X_1|^3] < +\infty\). Soit \(F_n\) la fonction de répartition de l’erreur renormalisée \(\frac{1}{\sigma \sqrt{n}} \sum_{k=1}^n X_k\). Alors

\[\begin{equation*} \forall x \in \mathbf{R}, \quad \forall n \ge 1, \quad |F_n(x) - F_{\mathcal{N}(0,1)}(x)| \le \frac{3 \rho}{\sigma^3 \sqrt{n}}. \end{equation*}\]

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

\[\begin{equation*} \forall r > 0, \quad \mathbf{P}\big(|X - \mathbf{E}[X]| \ge r \big) \le \frac{\operatorname{var}(X)}{r^2}. \end{equation*}\]

Appliquons cette proposition avec \(X = \frac{1}{n} \sum_{k=1}^n X_k\). Alors pour tout \(r > 0\)

\[\begin{equation*} \mathbf{P} \Big( \Big|\frac{1}{n} \sum_{k=1}^n X_k - m \Big| \ge r \Big) \le \frac{\frac{1}{n} \sum_{k=1}^n X_k}{r^2} = \frac{\sigma^2}{n r^2}, \end{equation*}\]

et

\[\begin{equation*} \forall n \ge 1, \, \forall r > 0, \quad \mathbf{P} \Big(m \in \Big[ \frac{1}{n}\sum_{k=1}^n X_k \pm r \Big] \Big) \ge 1 - \frac{\sigma^2}{n r^2}. \end{equation*}\]

Pour en déduire un I.C. de niveau \(1-\alpha\) il faut prendre \(r\) tel que \(\frac{\sigma^2}{n r^2} = \alpha\) c’est à dire \(r = \sqrt{1/\alpha} \frac{\sigma}{\sqrt{n}}\). Si on choisit \(\alpha = 0.05\) la constante \(\sqrt{1 / \alpha}\) qui apparait devant \(\frac{\sigma}{\sqrt{n}}\) est égale à \(4.47\). Pour le TCL on rappelle que la constante devant \(\frac{\sigma}{\sqrt{n}}\) vaut \(q_\alpha = 1.96\).

2.2.2. IC non asymptotique via Chernov

Si on suppose que \(X\) a des moments exponentiels, i.e. \(\exists \tau >0\), \(\mathbf{E}\big[ e^{\tau X} \big] < +\infty\) 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 \((X_n)_{n \ge 1}\) une suite de variables aléatoires réelles i.i.d. de loi \(\mu\) ayant des moments exponentiels et \(m = \mathbf{E}[X_1]\). Alors

\[\begin{align*} \forall a > m, \quad \mathbf{P} \Big(\frac{1}{n} \sum_{k=1}^n X_k > a \Big) &\le e^{-n h_\mu(a)}, \\ \forall a < m, \quad \mathbf{P} \Big(\frac{1}{n} \sum_{k=1}^n X_k < a \Big) &\le e^{-n h_\mu(a)}, \end{align*}\]

où \(h_\mu\) est la transformée de Cramer de la loi \(\mu\) définie par \(h_\mu(x) = \sup_{\theta \in \mathbf{R}} \big(\theta x - \psi(\theta)\big)\) et \(\psi(\theta) = \log \mathbf{E}\big[e^{\theta X}\big]\).

Appliquons cette proposition dans le cas Gaussien uniquement. On suppose donc que \(X_1 \sim \mathcal{N}(0, \sigma^2)\) (on prend \(m = 0\) pour simplifier). D’après cette proposition on a pour tout \(\varepsilon > 0\),

(2.7)\[\begin{equation} \mathbf{P} \Big( \Big|\frac{1}{n} \sum_{k=1}^n X_k\Big| > \varepsilon\Big) = \mathbf{P} \Big( \frac{1}{n} \sum_{k=1}^n X_k > \varepsilon \Big) + \mathbf{P} \Big(\frac{1}{n} \sum_{k=1}^n X_k < - \varepsilon \Big)\le 2 e^{-n h_\mu(\varepsilon)}. \end{equation}\]

La transformée de Cramer de la loi \(\mathcal{N}(0, \sigma^2)\) est donnée par \(h_\mu(x) = \frac{1}{2} \big(\frac{x}{\sigma}\big)^2\) donc

\[\begin{equation*} \mathbf{P} \Big( \Big| \frac{1}{n} \sum_{k=1}^n X_k \Big| \le \varepsilon \Big) \ge 1 - 2 e^{- \frac{n}{2} (\frac{\varepsilon}{\sigma})^2}. \end{equation*}\]

Pour construire un I.C. non asymptotique de niveau \(1-\alpha\) on choisit alors \(\varepsilon = \sqrt{-2 \log(\alpha/2)} \frac{\sigma}{\sqrt{n}}\). La constante pour \(\alpha = 0.05\) vaut alors \(\sqrt{-2 \log(0.025)} = 2.716\) (à comparer avec 1.96 et 4.47).