Simuler une binomiale

import numpy.random as rd
from math import comb

def bernoulli(p):
    return 1 if rd.random() < p else 0

def binomiale(n, p):
    # une binomiale EST une somme de n Bernoulli : on ecrit la definition
    return sum(bernoulli(p) for _ in range(n))

def comparer(n, p, tirages=100000):
    effectifs = [0] * (n + 1)
    for _ in range(tirages):
        effectifs[binomiale(n, p)] += 1
    for k in range(n + 1):
        empirique = effectifs[k] / tirages
        theorique = comb(n, k) * p**k * (1 - p)**(n - k)
        print(k, round(empirique, 4), round(theorique, 4))

Ce qu'on observe pour $n=10$, $p=0,3$ : les deux colonnes coïncident à environ $2\times 10^{-3}$ près, le mode étant en $k=3$ ($0,267$) et la moyenne empirique proche de $3,00$.

L'ordre de grandeur de l'écart. Chaque fréquence empirique est elle-même une variable de moyenne $P(X=k)$ et d'écart-type $\sqrt{P(1-P)/N}$. Pour $P\approx 0,27$ et $N=10^5$, cela donne $1,4\times 10^{-3}$ : les écarts constatés sont donc normaux, et non le signe d'une erreur. Un écart de $10^{-2}$, lui, mériterait qu'on relise le code.

Ce que la simulation prouve — et ne prouve pas. Elle conforte, elle chiffre, elle sert de garde-fou contre une erreur de calcul. Elle ne démontre rien : la loi des grands nombres garantit la convergence, pas l'exactitude d'un tirage fini.