Estimer une probabilité qu'on sait calculer

Le calcul exact. Conditionnons par la valeur de $X$. Pour $x\in ]0,1]$, l'événement $[XY⩽\frac{1}{2}]$ s'écrit $[Y⩽\frac{1}{2x}]$, dont la probabilité vaut $\min \left(1,\frac{1}{2x}\right)$ puisque $Y$ est uniforme sur $[0,1]$ et indépendante de $X$. Cette quantité vaut $1$ tant que $x⩽\frac{1}{2}$. En intégrant sur la densité de $X$, constante égale à $1$ sur $[0,1]$ : $$p={\int }_0^{1/2}1\mathrm{d}x+{\int }_{1/2}^1\frac{1}{2x}\mathrm{d}x=\frac{1}{2}+\frac{1}{2}[\ln x{]}_{1/2}^1=\frac{1}{2}+\frac{\ln 2}{2}\approx 0,8466.$$

L'estimateur de Monte-Carlo. On tire $n$ couples et l'on note $F_n$ la proportion de ceux qui vérifient $XY⩽\frac{1}{2}$. Chaque indicatrice suit une loi de Bernoulli de paramètre $p$, donc $E(F_n)=p$ : l'estimateur est sans biais, et convergent par la loi faible des grands nombres.

La demi-largeur. La variance d'une indicatrice vaut $p(1-p)$, ici $$0,8466\times 0,1534\approx 0,1299,\sigma \approx 0,3604.$$ Pour $n=10^6$, le théorème limite central donne $$\frac{1,96\sigma }{\sqrt{n}}=\frac{1,96\times 0,3604}{1000}\approx 7,1\times 10^{-4}.$$ On connaît donc $p$ à environ $0,0007$ près, soit trois décimales sûres.

Le code correspondant.

import numpy as np
rng = np.random.default_rng(3)
n = 10**6
X = rng.random(n)
Y = rng.random(n)
F = float(np.mean(X*Y <= 0.5))
print(round(F, 6), round(1.96*np.sqrt(F*(1 - F)/n), 6))

Pourquoi commencer par un cas calculable. On ne peut faire confiance à un code que si on l'a d'abord confronté à une valeur connue : ici $\frac{1}{2}+\frac{\ln 2}{2}$ sert de témoin.