Estimer une aire par Monte-Carlo

La probabilité visée. Le couple $(U,V)$ est uniforme sur un carré d'aire $1$, donc la probabilité de tomber dans le quart de disque vaut l'aire de celui-ci : $p=\frac{\pi }{4}\approx 0,7854$.

L'estimateur. $n{F}_n$ compte les succès de $n$ épreuves indépendantes de probabilité $p$, donc $n{F}_n\hookrightarrow \mathcal{B}(n,p)$ et $E(4F_n)=4\times \frac{\pi }{4}=\pi$ : l'estimateur est sans biais, pour tout $n$.

Sa variance. $V(F_n)=\frac{p(1-p)}{n}$, donc $V(4F_n)=\frac{16p(1-p)}{n}\to 0$ : avec l'absence de biais, la condition suffisante donne la convergence.

La demi-largeur. $p(1-p)=0,7854\times 0,2146=0,1686$, et pour $n=10^6$ :

$$4\times 1,96\sqrt{\frac{0,1686}{10^6}}=4\times 1,96\times 4,11\times 10^{-4}\approx 3,2\times 10^{-3}.$$

Un million de tirages donnent donc $\pi$ à $\pm 0,0032$ : deux décimales sûres, pas trois.

Le code.

import numpy as np
import numpy.random as rd
n = 10**6
u, v = rd.random(n), rd.random(n)
f = np.mean(u**2 + v**2 <= 1)
print(4*f, 4*1.96*np.sqrt(f*(1-f)/n))

Le point à retenir. Monte-Carlo converge en $\frac{1}{\sqrt{n}}$ et rien ne l'accélère : une décimale de plus coûte cent fois plus de tirages. C'est pourquoi on ne calcule jamais $\pi$ ainsi — mais c'est aussi pourquoi la méthode reste imbattable pour ce qu'aucune formule ne donne.