Estimer $\pi$ par Monte-Carlo

import numpy.random as rd

def estime_pi(n):
    dedans = 0
    for _ in range(n):
        x, y = rd.random(), rd.random()
        if x * x + y * y <= 1:
            dedans += 1
    return 4 * dedans / n

Le principe. Le quart de disque a pour aire $\frac{\pi }{4}$ et le carré $1$. La proportion de points tombés dedans estime donc $\frac{\pi }{4}$, qu'on multiplie par $4$.

Combien de tirages. L'estimateur a pour écart-type $\frac{4\sqrt{p(1-p)}}{\sqrt{n}}$ avec $p=\frac{\pi }{4}$, soit environ $\frac{1,64}{\sqrt{n}}$. Pour espérer deux décimales, il faut cet écart-type sous $5\times 10^{-3}$, donc

$$\sqrt{n}>\frac{1,64}{0,005}\approx 328,n>10^5\text{environ}.$$

Et encore : à $n=10^5$ on obtient typiquement $3,14$ ou $3,15$ — les décimales ne sont pas stables, elles fluctuent d'une exécution à l'autre.

Pourquoi c'est une mauvaise méthode pour $\pi$. L'erreur décroît en $\frac{1}{\sqrt{n}}$ : gagner une décimale demande de multiplier le nombre de tirages par cent. Pour dix décimales il faudrait $10^{20}$ points. La méthode des trapèzes du chapitre 8, avec son erreur en $n^{-2}$, fait mieux avec mille points ; et la série $\pi =4\sum \frac{(-1{)}^k}{2k+1}$, ou mieux les formules de Machin, donnent des milliers de décimales en quelques secondes.

Alors pourquoi Monte-Carlo existe. Parce que son erreur en $\frac{1}{\sqrt{n}}$ ne dépend pas de la dimension. En dimension $20$, une grille de trapèzes est hors de portée ; Monte-Carlo, lui, marche exactement pareil. C'est là qu'il est irremplaçable — et $\pi$, en dimension deux, est justement le cas où il est le plus mauvais.