Estimer une aire en tirant des points

Le principe. Un point $(x;y)$ tiré au hasard dans le carré unité tombe dans une région avec une probabilité égale à l'aire de cette région, le carré étant d'aire $1$. Or la région située sous la courbe de $x\mapsto x^2$ et au-dessus de l'axe a pour aire ${\int }_0^1{x}^2\mathrm{d}x$. La fréquence des points vérifiant $y\le x^2$ estime donc cette intégrale.

Le script.

import numpy.random as rd

n = 10000
compte = 0
for i in range(n):
    x = rd.random()
    y = rd.random()
    if y <= x**2:
        compte = compte + 1
print(compte/n)   # 0.3327 par exemple

La valeur exacte. Une primitive de $x^2$ est $\frac{x^3}{3}$, donc $${\int }_0^1{x}^2\mathrm{d}x={\left[\frac{x^3}{3}\right]}_0^1=\frac{1}{3}\approx 0,3333.$$

La comparaison. La fréquence obtenue, environ $0,333$, en est très proche ; elle varie de quelques millièmes d'une exécution à l'autre.

Le point délicat. Les deux appels à rd.random() doivent être séparés : écrire x = rd.random() puis y = x tirerait toujours des points de la diagonale et donnerait un résultat sans rapport. Chaque appel produit un nouveau nombre, indépendant du précédent — c'est ce qui fait des deux coordonnées un point réellement au hasard dans le carré.