Une fréquence, et la probabilité exacte

Ce qu'on montre. Que la simulation et le calcul exact se contrôlent l'un l'autre, et que la loi binomiale se calcule par récurrence, sans coefficient binomial explicite.

La récurrence. Pour $X\hookrightarrow \mathcal{B}(n,p)$ et $0\le k<n$, $$\frac{P(X=k+1)}{P(X=k)}=\frac{\left(\genfrac{}{}{0ex}{}{n}{k+1}\right)p^{k+1}(1-p{)}^{n-k-1}}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k}}=\frac{n-k}{k+1}\cdot \frac{p}{1-p},$$ en simplifiant les factorielles et les puissances. On part donc de $P(X=0)=(1-p{)}^n$ et l'on avance de proche en proche.

Le programme.

import numpy as np
import numpy.random as rd

# --- par simulation
X = rd.binomial(10, 0.3, 10 ** 6)
print(np.mean(X >= 3))          # ~0.6172

# --- exactement
def proba_au_moins(n, p, m):
    """P(X >= m) pour X binomiale(n, p), par recurrence sur la loi."""
    q = (1 - p) ** n            # P(X = 0)
    cumul = q
    for k in range(m - 1):      # ajoute P(X=1), ..., P(X=m-1)
        q = q * (n - k) / (k + 1) * p / (1 - p)
        cumul = cumul + q
    return 1 - cumul

print(proba_au_moins(10, 0.3, 3))   # 0.617219...

Le calcul à la main. $$P(X=0)=0,7^{10}\approx 0,02825,P(X=1)=10\times 0,3\times 0,7^9\approx 0,12106,$$ $$P(X=2)=45\times 0,09\times 0,7^8\approx 0,23347.$$ Leur somme vaut $\approx 0,38278$, donc $$P(X\ge 3)=1-0,38278\approx 0,61722.$$

La comparaison. La simulation sur $10^6$ tirages affiche une valeur voisine de $0,6172$. L'écart-type de cette fréquence vaut $\sqrt{\frac{0,617\times 0,383}{10^6}}\approx 4,9\cdot 10^{-4}$ : les deux résultats doivent coïncider sur trois décimales, ce qu'on observe.

Point délicat. La boucle parcourt range(m - 1), soit $k=0,\dots ,m-2$, ce qui ajoute $P(X=1),\dots ,P(X=m-1)$ à la valeur initiale $P(X=0)$ : on cumule bien $P(X\le m-1)$. Écrire range(m) ajouterait un terme de trop et renverrait $P(X\ge m+1)$ — erreur silencieuse, que seule la comparaison avec la simulation révélerait.