Comprenez la loi des grands nombres

Vous avez écrit un programme non déterministe : un algorithme randomisé dont le temps d’exécution change à chaque lancement. Une première exécution prend 2,0 secondes. Que pouvez-vous en conclure sur son temps « typique » ? Pas grand-chose. Vous le relancez trois fois (1,4 s, 1,9 s, 1,7 s) et vous gardez la moyenne des quatre mesures, 1,75 s. Instinctivement, vous faites plus confiance à cette moyenne. Mais pourquoi, exactement ? Et combien d’exécutions faut-il pour une précision donnée ?

La réponse tient en un théorème, la loi des grands nombres, qui justifie à la fois la vision fréquentiste du chapitre 1 et toutes les simulations de ce cours. Le fil rouge sera la méthode de Monte-Carlo, qui estime une probabilité ou une intégrale par simulation.

La moyenne d’un échantillon

Répéter une expérience $n$ fois dans les mêmes conditions, c’est observer $n$ variables $X_1, \dots, X_n$ indépendantes et identiquement distribuées (on écrit i.i.d.) : même loi que $X$, et aucune influence de l’une sur l’autre. On note $\mu = E[X]$ et $\sigma^2 = \mathrm{Var}(X)$, et l’on forme la moyenne empirique :

$$\bar{X}_n = \frac{1}{n} \sum_{i=1}^{n} X_i$$

$\bar{X}_n$ est elle-même une variable aléatoire : relancez vos $n$ exécutions, vous obtiendrez une autre moyenne. Les chapitres 4 et 6 donnent immédiatement ses deux premiers moments.

• Par linéarité de l’espérance : $E[\bar{X}_n] = \frac{1}{n} \times n\mu = \mu$. En moyenne, la moyenne empirique tombe juste.
• Par indépendance, les variances s’additionnent : $\mathrm{Var}(\bar{X}_n) = \frac{1}{n^2} \times n\sigma^2 = \frac{\sigma^2}{n}$.

La seconde formule est la clé de tout ce qui suit : la moyenne varie moins qu’une mesure seule, et sa variance tend vers 0 quand $n$ grandit. Son écart-type, $\sigma / \sqrt{n}$, ne décroît cependant qu’en $1/\sqrt{n}$ : pour diviser l’incertitude par 10, il faut 100 fois plus de mesures.

Toutes ces lois sont centrées sur μ = 3,5, mais plus n grandit, plus elles se resserrent : quadrupler le nombre de dés divise l'écart-type par deux.

Borner les écarts : l’inégalité de Bienaymé-Tchebychev

Une variance faible signifie que la variable s’écarte rarement de son espérance. L’inégalité de Bienaymé-Tchebychev quantifie ce « rarement », pour n’importe quelle loi de variance finie :

$$P\big(\lvert X - E[X] \rvert \geq \varepsilon\big) \leq \frac{\mathrm{Var}(X)}{\varepsilon^2} \quad \text{pour tout } \varepsilon \gt 0$$

• $\lvert X - E[X] \rvert$ est l’écart (en valeur absolue) entre $X$ et son espérance ;
• $\varepsilon$ est le seuil d’écart qui vous inquiète.

En posant $\varepsilon = k\sigma$, elle devient : la probabilité de s’écarter de plus de $k$ écarts-types est au plus $1/k^2$. Par exemple, au plus 25 % de chances de s’écarter de plus de $2\sigma$, quelle que soit la loi.

Exemple

Pour une normale, la probabilité de s’écarter de plus de $2\sigma$ vaut 0,046 (chapitre 5). Tchebychev garantit seulement « au plus 0,25 ». La borne est très pessimiste, parce qu’elle doit valoir pour toutes les lois, même les plus étalées. Son intérêt n’est pas sa précision, mais sa généralité.

La borne de Tchebychev reste toujours au-dessus de la vraie probabilité ; pour une normale, elle la surestime d'un facteur 5 dès k = 2.

L’inégalité se démontre en deux lignes à partir d’une inégalité encore plus simple, celle de Markov : pour $Y \geq 0$ et $a \gt 0$, $P(Y \geq a) \leq E[Y]/a$. (En effet, $Y \geq a \, \mathbb{1}_{Y \geq a}$, dont on prend l’espérance.) On l’applique à $Y = (X - E[X])^2$ et $a = \varepsilon^2$.

La loi des grands nombres

Appliquons Tchebychev à la moyenne empirique, de variance $\sigma^2/n$ :

$$P\big(\lvert \bar{X}_n - \mu \rvert \geq \varepsilon\big) \leq \frac{\sigma^2}{n \varepsilon^2} \xrightarrow[n \to \infty]{} 0$$

C’est la loi faible des grands nombres : quel que soit le seuil $\varepsilon$, la probabilité que la moyenne empirique s’écarte de $\mu$ de plus de $\varepsilon$ tend vers 0. On dit que $\bar{X}_n$ converge en probabilité vers $\mu$. (Une version plus forte, la loi forte des grands nombres, affirme même que la suite des moyennes converge vers $\mu$ avec probabilité 1 ; elle ne demande que l’existence de l’espérance.)

Chaque courbe est une suite de lancers. Au début, la moyenne erre entre 1 et 6 ; ensuite, toutes les trajectoires sont prises dans l'entonnoir, qui se referme en 1/√n.

À retenir

Pour des $X_i$ i.i.d. d’espérance $\mu$ finie, la moyenne empirique $\bar{X}_n$ converge vers $\mu$. Elle est sans biais ($E[\bar{X}_n] = \mu$) et sa variance vaut $\sigma^2/n$ : l’incertitude décroît en $\sigma/\sqrt{n}$.

La loi des grands nombres justifie enfin la vision fréquentiste du chapitre 1. Pour un événement $A$ de probabilité $p$, posons $Y_i = \mathbb{1}_{A}$ lors de la $i$-ème répétition (l’indicatrice : 1 si $A$ se réalise, 0 sinon). Les $Y_i$ suivent la loi de Bernoulli de paramètre $p$, et leur moyenne est exactement la fréquence de $A$. Elle converge donc vers $E[Y_i] = p$. Probabilité et fréquence à long terme coïncident : c’est un théorème, pas une définition.

Quand la loi des grands nombres échoue : la loi de Cauchy

L’hypothèse « espérance finie » n’est pas une formalité. La loi de Cauchy, de densité $f(x) = \frac{1}{\pi(1 + x^2)}$, a une allure de cloche symétrique autour de 0, un peu comme une normale. Mais ses queues décroissent si lentement que $\int \lvert x \rvert f(x) \, dx$ diverge : elle n’a pas d’espérance (ni de variance). Conséquence spectaculaire : la moyenne de $n$ variables de Cauchy indépendantes suit encore une loi de Cauchy, exactement aussi dispersée qu’une seule observation. Moyenner ne sert à rien.

set.seed(8)
n <- 10000
x_exp    <- rexp(n, rate = 1)            # espérance 1
x_cauchy <- rcauchy(n)                   # pas d'espérance

par(mfrow = c(1, 2))
plot(cumsum(x_exp) / (1:n), type = "l", ylab = "moyenne cumulée", main = "Exponentielle")
abline(h = 1, lty = 2)                   # la moyenne se stabilise autour de 1
plot(cumsum(x_cauchy) / (1:n), type = "l", ylab = "moyenne cumulée", main = "Cauchy")
# la moyenne fait des sauts brutaux à chaque valeur extrême, et ne se stabilise jamais

Pour l'exponentielle, la moyenne se stabilise vite. Pour la loi de Cauchy, une seule valeur extrême suffit à faire sauter la moyenne, même après 8 000 tirages : elle ne converge pas.

La loi de Cauchy n’est pas une curiosité : on l’obtient comme rapport de deux normales centrées indépendantes, ou comme la position où frappe un rayon lumineux émis dans une direction aléatoire. Plus généralement, les données à queues lourdes (revenus, tailles de fichiers, temps de réponse extrêmes) rendent la moyenne très instable : pensez alors à la médiane.

Piège

La loi des grands nombres suppose une espérance finie. Si vos données contiennent des valeurs extrêmes qui dominent régulièrement la somme, la moyenne empirique peut ne jamais se stabiliser. Tracez toujours la moyenne cumulée : si elle fait des sauts qui ne s’amortissent pas, méfiez-vous.

Le fil rouge : la méthode de Monte-Carlo

La loi des grands nombres a une application immédiate : pour calculer une espérance $E[g(X)]$ (et donc une probabilité, qui est l’espérance d’une indicatrice), il suffit de simuler beaucoup de tirages et d’en faire la moyenne. C’est la méthode de Monte-Carlo, du nom du casino. Elle sert quand le calcul exact est impossible ou pénible : intégrales en grande dimension, systèmes complexes, lois de sommes.

Estimer π avec des fléchettes

Tirez un point $(X, Y)$ uniformément dans le carré $[-1\,;\,1]^2$. La probabilité qu’il tombe dans le disque unité est le rapport des aires, $\pi/4$. Donc $\pi = 4 \, P(X^2 + Y^2 \leq 1)$, que l’on estime par une fréquence :

set.seed(1)
n <- 1e6
x <- runif(n, -1, 1); y <- runif(n, -1, 1)
dedans <- (x^2 + y^2 <= 1)               # indicatrice « dans le disque »
4 * mean(dedans)                         # estimation de pi
4 * sd(dedans) / sqrt(n)                 # ordre de grandeur de l'erreur : environ 0.0016

L’erreur typique est l’écart-type de l’estimateur, $4\sqrt{p(1-p)/n}$ avec $p = \pi/4$ : environ $0{,}016$ pour $n = 10^4$ tirages, et $0{,}0016$ pour $n = 10^6$. Cent fois plus de tirages pour un chiffre de plus : c’est la loi en $1/\sqrt{n}$, le prix à payer de Monte-Carlo.

La fréquence des points dans le disque estime π/4. L'estimation reste dans une bande de largeur 4√(p(1 − p)/n) : il faut passer de 10 000 à 1 000 000 de points pour diviser l'erreur par 10.

Estimer une probabilité sans la calculer

Au chapitre 6, la durée totale pour aller au campus était $U + T$, avec $U \sim \mathcal{U}[0, 10]$ et $T \sim \mathcal{N}(30, 4^2)$ indépendantes, et nous avions dû calculer une convolution. Avec Monte-Carlo, deux lignes suffisent :

set.seed(6)
n <- 1e6
duree <- runif(n, 0, 10) + rnorm(n, 30, 4)
mean(duree > 45)                         # P(être en retard avec 45 min de marge), environ 0.020
quantile(duree, 0.95)                    # environ 43.1 minutes

Combien de tirages ?

Tchebychev répond à la question « combien de simulations pour une précision donnée ? ». Pour estimer une probabilité $p$ à $\pm 0{,}01$ près avec une probabilité d’erreur d’au plus 5 %, on veut $\frac{p(1-p)}{n \times 0{,}01^2} \leq 0{,}05$. Comme $p(1 - p) \leq 1/4$ pour tout $p$, il suffit que :

$$n \geq \frac{1}{4 \times 0{,}01^2 \times 0{,}05} = 50\,000$$

Ce nombre est une garantie, valable pour toute valeur de $p$. Le théorème central limite, au chapitre suivant, montrera qu’environ 10 000 tirages suffisent en réalité : Tchebychev est prudent.

En entretien

« La loi des grands nombres dit que la moyenne de tirages i.i.d. converge vers l’espérance, pourvu qu’elle existe. C’est ce qui justifie les simulations de Monte-Carlo, dont l’erreur décroît en $1/\sqrt{n}$ : pour gagner un chiffre de précision, il faut cent fois plus de tirages. »

En résumé

• Pour $n$ variables i.i.d., la moyenne empirique $\bar{X}_n$ a pour espérance $\mu$ et pour variance $\sigma^2/n$.
• L’inégalité de Bienaymé-Tchebychev $P(\lvert X - E[X] \rvert \geq \varepsilon) \leq \mathrm{Var}(X)/\varepsilon^2$ vaut pour toute loi de variance finie, au prix d’une borne pessimiste.
• La loi des grands nombres : $\bar{X}_n$ converge vers $\mu$ ; appliquée aux indicatrices, elle montre que la fréquence converge vers la probabilité.
• Elle exige une espérance finie : la loi de Cauchy en est le contre-exemple, sa moyenne empirique ne se stabilise jamais.
• La méthode de Monte-Carlo estime une espérance, une probabilité ou une intégrale par une moyenne de simulations, avec une erreur en $1/\sqrt{n}$.

Exercices

Exercice · Borner la moyenne de 100 dés

facile

On lance 100 dés équilibrés et on calcule la moyenne $\bar{X}_{100}$ des faces. Majorez la probabilité que cette moyenne s’écarte de 3,5 d’au moins 0,5 avec l’inégalité de Bienaymé-Tchebychev.

Voir la solution

Pour un dé, $\sigma^2 = 35/12$ (chapitre 4). Donc $\mathrm{Var}(\bar{X}_{100}) = \frac{35/12}{100} \approx 0{,}029$, et :

$$P\big(\lvert \bar{X}_{100} - 3{,}5 \rvert \geq 0{,}5\big) \leq \frac{0{,}029}{0{,}5^2} \approx 0{,}117$$

Au plus 12 % environ. La vraie probabilité est bien plus faible (de l’ordre de $0{,}3\,\%$, comme le montrera le théorème central limite), ce qui illustre le caractère prudent de la borne.

En code

Python

import numpy as np

var_de = 35 / 12
var_moyenne = var_de / 100
print(var_moyenne / 0.5**2)                 # 0.1167  borne de Tchebychev

# Vraie probabilité, estimée sur 200 000 séries de 100 dés
rng = np.random.default_rng(0)
moyennes = rng.integers(1, 7, size=(200_000, 100)).mean(axis=1)
print(np.mean(np.abs(moyennes - 3.5) >= 0.5))   # 0.0036, bien sous la borne

R

var_de <- 35 / 12
var_moyenne <- var_de / 100
var_moyenne / 0.5^2                         # 0.1167  borne de Tchebychev

# Vraie probabilité, estimée sur 200 000 séries de 100 dés
set.seed(0)
des <- matrix(sample(1:6, 100 * 2e5, replace = TRUE), ncol = 100)
moyennes <- rowMeans(des)
mean(abs(moyennes - 3.5) >= 0.5)            # 0.0037, bien sous la borne

Exercice · Une intégrale par Monte-Carlo

moyen

On veut calculer $I = \int_0^1 e^{-x^2} \, dx$, qui n’a pas de primitive simple.

1. Montrez que $I = E[g(U)]$ pour une fonction $g$ et une variable $U$ que vous préciserez.
2. Écrivez le code R qui estime $I$ avec $10^5$ tirages.
3. On admet que $\mathrm{Var}(g(U)) \approx 0{,}040$. Combien de tirages garantissent, par Tchebychev, une erreur inférieure à $0{,}001$ avec probabilité au moins 95 % ?

Voir la solution

1. Si $U \sim \mathcal{U}[0, 1]$, sa densité vaut 1 sur $[0, 1]$, donc par le théorème de transfert $E[e^{-U^2}] = \int_0^1 e^{-x^2} \times 1 \, dx = I$. On prend $g(x) = e^{-x^2}$.
2. Le code est donné plus bas, dans « En code » : avec $10^5$ tirages, on obtient environ $0{,}747$, pour une valeur exacte de $0{,}7468$.
3. On veut $\frac{0{,}040}{n \times 0{,}001^2} \leq 0{,}05$, soit $n \geq \frac{0{,}040}{0{,}05 \times 10^{-6}} = 800\,000$ tirages. C’est beaucoup pour une intégrale en dimension 1 (une méthode des trapèzes ferait mieux), mais le nombre de tirages de Monte-Carlo ne dépend pas de la dimension, ce qui la rend imbattable pour les intégrales en dimension 10 ou 100.

En code

Python

import numpy as np
from scipy.integrate import quad

rng = np.random.default_rng(9)
u = rng.random(100_000)                     # U uniforme sur [0, 1]
g_u = np.exp(-u**2)
print(g_u.mean())                           # 0.7478
print(quad(lambda x: np.exp(-x**2), 0, 1)[0])   # 0.7468 (valeur de référence)

# 3. Tirages garantis par Tchebychev : Var / (n eps^2) <= 0.05
print(g_u.var())                            # 0.0404
print(0.040 / (0.05 * 0.001**2))            # 800000

R

set.seed(9)
u <- runif(1e5)                             # U uniforme sur [0, 1]
mean(exp(-u^2))                             # 0.7471
integrate(function(x) exp(-x^2), 0, 1)      # 0.7468 (valeur de référence)

# 3. Tirages garantis par Tchebychev : Var / (n eps^2) <= 0.05
var(exp(-u^2))                              # 0.0403
0.040 / (0.05 * 0.001^2)                    # 8e+05

Exercice · Démontrer Markov et Tchebychev

difficile

Soit $Y$ une variable aléatoire positive, d’espérance finie, et $a \gt 0$.

1. Montrez que $Y \geq a \, \mathbb{1}_{Y \geq a}$, puis que $P(Y \geq a) \leq E[Y]/a$ (inégalité de Markov).
2. Déduisez-en l’inégalité de Bienaymé-Tchebychev.
3. Un serveur a un temps de réponse moyen de 200 ms. Que peut-on dire, sans autre information, de la proportion de requêtes qui prennent plus d’une seconde ? Et si l’on sait en plus que l’écart-type vaut 100 ms ?

Voir la solution

1. Si $Y \geq a$, l’indicatrice vaut 1 et l’inégalité devient $Y \geq a$, vraie. Sinon, elle vaut 0 et l’inégalité devient $Y \geq 0$, vraie car $Y$ est positive. En prenant l’espérance (qui préserve les inégalités) : $E[Y] \geq a \, E[\mathbb{1}_{Y \geq a}] = a \, P(Y \geq a)$, d’où le résultat.
2. On applique Markov à $Y = (X - E[X])^2 \geq 0$ et $a = \varepsilon^2$ : $P\big((X - E[X])^2 \geq \varepsilon^2\big) \leq E[(X - E[X])^2]/\varepsilon^2 = \mathrm{Var}(X)/\varepsilon^2$. Or $(X - E[X])^2 \geq \varepsilon^2$ équivaut à $\lvert X - E[X] \rvert \geq \varepsilon$.
3. Par Markov (le temps est positif) : $P(T \geq 1\,000) \leq 200/1\,000 = 0{,}2$. Au plus 20 % des requêtes. Avec l’écart-type, Tchebychev donne mieux : $T \geq 1\,000$ implique $\lvert T - 200 \rvert \geq 800$, donc $P(T \geq 1\,000) \leq 100^2/800^2 \approx 0{,}016$, soit au plus 1,6 %. Plus on connaît d’information sur la loi, plus la borne est fine.

En code

Python

from scipy import stats

# Bornes sans connaître la loi
print(200 / 1000)                 # 0.2      Markov : P(T >= 1000) <= E[T] / 1000
print(100**2 / 800**2)            # 0.015625 Tchebychev : écart d'au moins 800 ms

# Exemple de loi compatible : gamma de moyenne 200 et d'écart-type 100
t = stats.gamma(a=4, scale=50)
print(t.mean(), t.std())          # 200.0 100.0
print(t.sf(1000))                 # 3.2e-06 : les bornes sont très larges

R

# Bornes sans connaître la loi
200 / 1000                        # 0.2      Markov : P(T >= 1000) <= E[T] / 1000
100^2 / 800^2                     # 0.015625 Tchebychev : écart d'au moins 800 ms

# Exemple de loi compatible : gamma de moyenne 200 et d'écart-type 100
c(4 * 50, sqrt(4 * 50^2))                 # 200 100
1 - pgamma(1000, shape = 4, scale = 50)   # 3.2e-06 : les bornes sont très larges