Estimer $\pi$ en tirant des points au hasard

Ce qu'on montre. Que l'estimation d'une probabilité par une fréquence est fondée par la loi faible des grands nombres, et que Bienaymé-Tchebychev en chiffre le coût.

Le fondement. Notons $A$ l'événement $[X^2+Y^2\le 1]$. Le point $(X,Y)$ étant uniforme sur le carré, la probabilité de $A$ est le rapport de l'aire du quart de disque à celle du carré : $$P(A)=\frac{\pi /4}{1}=\frac{\pi }{4}.$$ Posons $Z_i=1$ si le $i$-ième point tombe dans le quart de disque et $Z_i=0$ sinon : les $Z_i$ sont indépendantes, de même loi de Bernoulli de paramètre $p=\frac{\pi }{4}$. La loi faible des grands nombres affirme que $\overline{Z_N}\to p$ au sens du cours, donc $4\overline{Z_N}$ estime $\pi$.

Le programme.

import numpy as np
import numpy.random as rd

def estime_pi(N):
    X = rd.random(N)                 # N uniformes sur [0,1[
    Y = rd.random(N)
    dedans = (X ** 2 + Y ** 2) <= 1  # tableau de booleens
    return 4 * np.mean(dedans)

for N in [10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6]:
    print(N, estime_pi(N))

La moyenne d'un tableau de booléens est exactement la fréquence des True : c'est l'estimation cherchée.

La précision. Posons ${\hat{\pi }}_N=4\overline{Z_N}$. Alors $E({\hat{\pi }}_N)=4p=\pi$, et $$V\left({\hat{\pi }}_N\right)=16V\left(\overline{Z_N}\right)=\frac{16p(1-p)}{N},$$ la variance d'une moyenne de $N$ variables indépendantes de même loi valant $\frac{V(Z_1)}{N}$. Numériquement $p=\frac{\pi }{4}\approx 0,7854$, donc $p(1-p)\approx 0,1685$ et $$V\left({\hat{\pi }}_N\right)\approx \frac{2,696}{N}.$$

Le $N$ demandé. Bienaymé-Tchebychev, appliquée à ${\hat{\pi }}_N$ avec $\epsilon =0,01$ : $$P(|{\hat{\pi }}_N-\pi |\ge 0,01)\le \frac{2,696}{N\times 10^{-4}}.$$ Pour que ce majorant soit inférieur à $0,05$, il suffit que $$N\ge \frac{2,696}{10^{-4}\times 0,05}\approx 5,4\cdot 10^5.$$ Environ 540 000 tirages.

Commentaire. La précision progresse comme $\frac{1}{\sqrt{N}}$ : gagner une décimale demande de multiplier $N$ par cent — c'est pourquoi on ne calcule pas $\pi$ ainsi. La force de la méthode est ailleurs : le nombre de tirages ne dépend que de $\epsilon$, jamais de la complexité de l'événement.