Simuler des lois de probabilité en Python et au tableur

Simuler une expérience aléatoire, c'est la reproduire un grand nombre de fois avec des nombres tirés au hasard par l'ordinateur, puis mesurer la fréquence de l'événement qui nous intéresse : par la loi des grands nombres, elle se rapproche de sa probabilité. Le programme du BTS insiste sur cette démarche pour chaque loi : simuler un schéma de Bernoulli, une loi uniforme, et obtenir les lois exponentielle et normale à partir de la loi uniforme. Ce chapitre le fait en Python, avec la seule bibliothèque standard, puis montre les calculs exacts et leur équivalent au tableur.

Pourquoi simuler

Une simulation sert à trois choses :

• vérifier un calcul : si la formule donne 0,0115 et que la simulation donne 0,0498, l'un des deux est faux ;
• comprendre une loi : voir l'histogramme d'une somme d'uniformes prendre la forme d'une cloche convainc mieux qu'un théorème ;
• estimer ce qu'on ne sait pas calculer : la durée moyenne d'indisponibilité d'une infrastructure compliquée, le temps d'attente dans une file avec plusieurs techniciens.

Tous les programmes de ce chapitre s'exécutent avec Python 3.12 ou plus récent, sans rien installer : lance-les dans un terminal ou dans un éditeur comme VS Code. L'application Python des calculatrices n'en fait tourner qu'une partie : son module random est réduit (ni choices, ni expovariate, ni gauss), elle n'a ni statistics ni math.comb, et elle est bien plus lente qu'un ordinateur. Sur calculatrice, garde les simulations à base de random.random() et réduis n à 1 000.

Le module random

Estimer une probabilité et juger la précision

Le schéma est toujours le même : répéter n fois l'expérience, compter les réalisations de l'événement, diviser par n.

import random

def frequence(experience, n):
    """Répète n fois experience(), qui renvoie True ou False, et renvoie la fréquence des True."""
    return sum(1 for _ in range(n) if experience()) / n

def aucun_paquet_perdu():
    return all(random.random() >= 0.05 for _ in range(10))   # 10 paquets, perte 5 %

print(frequence(aucun_paquet_perdu, 100_000))   # environ 0.599 (exact : 0,95¹⁰ ≈ 0,5987)

Quelle confiance accorder à cette fréquence ? Une règle simple, vue au lycée : avec n répétitions, la fréquence observée tombe dans 95 % des cas à moins de 1/√n de la probabilité. Avec 100 000 répétitions, 1/√100 000 ≈ 0,003 : on peut faire confiance aux deux premières décimales. Pour gagner une décimale, il faut 100 fois plus de répétitions.

Simuler un schéma de Bernoulli et une loi binomiale

Une épreuve de Bernoulli de paramètre p se simule par random.random() < p : le nombre uniforme tombe sous p avec la probabilité p. Le nombre de succès en n épreuves s'obtient en comptant :

import random
from math import comb

def binomiale(n, p):
    """Un tirage de la loi B(n ; p) : le nombre de succès en n épreuves de Bernoulli."""
    return sum(1 for _ in range(n) if random.random() < p)

def proba_binomiale(n, p, k):
    """P(X = k) exacte pour X de loi B(n ; p)."""
    return comb(n, k) * p**k * (1 - p) ** (n - k)

random.seed(2026)
N = 100_000
tirages = [binomiale(30, 0.2) for _ in range(N)]

print("moyenne :", sum(tirages) / N)              # proche de np = 6
for k in range(0, 13):
    freq = tirages.count(k) / N
    barre = "#" * round(freq * 200)
    print(f"{k:2d}  {freq:.4f}  {proba_binomiale(30, 0.2, k):.4f}  {barre}")

Le programme affiche, pour chaque valeur k, la fréquence simulée, la probabilité exacte et une barre proportionnelle : un diagramme en bâtons couché, qui reproduit celui du chapitre 5. La moyenne des tirages est proche de 6, ce qui conforte la formule E(X) = np que le programme admet.

Simuler une loi discrète quelconque

Pour une loi donnée par un tableau, random.choices tire une valeur en respectant les probabilités :

import random

valeurs = [0, 1, 2, 3, 4]                 # tickets urgents reçus en une heure
probas = [0.30, 0.35, 0.20, 0.10, 0.05]

heures = random.choices(valeurs, weights=probas, k=100_000)
print(sum(heures) / len(heures))          # proche de E(X) = 1,25

Sur une journée de 8 heures, sum(random.choices(valeurs, weights=probas, k=8)) simule le nombre de tickets urgents de la journée : en répétant ce tirage, tu estimes par exemple la probabilité qu'une journée en compte plus de 15, sans calculer la loi de la somme.

Simuler la loi exponentielle à partir de la loi uniforme

Si U suit la loi uniforme sur [0 ; 1[, alors T = −ln(1 − U) / λ suit la loi exponentielle de paramètre λ. La preuve tient en une ligne : pour t ≥ 0,

P(T ≤ t) = P(−ln(1 − U) ≤ λt) = P(1 − U ≥ e^(−λt)) = P(U ≤ 1 − e^(−λt)) = 1 − e^(−λt)

car pour une loi uniforme sur [0 ; 1[, P(U ≤ x) = x. On retrouve la fonction de répartition de la loi exponentielle. Cette méthode, dite de la fonction de répartition inverse, transforme un nombre uniforme en n'importe quelle loi dont on sait inverser la fonction de répartition.

import random
from math import log, exp

def exponentielle(lam):
    return -log(1 - random.random()) / lam   # 1 - random() est dans ]0 ; 1], log ne plante jamais

lam = 1 / 100_000                            # alimentation de MTBF 100 000 h
N = 100_000
durees = [exponentielle(lam) for _ in range(N)]
print(sum(durees) / N)                                 # proche de 100 000
print(sum(1 for d in durees if d <= 8760) / N)         # proche de 1 - exp(-0.0876) ≈ 0.0839
print(1 - exp(-lam * 8760))

Simuler une loi de Poisson avec des attentes exponentielles

Le programme introduit la loi de Poisson comme le nombre d'événements observés pendant une durée donnée quand les attentes entre événements suivent une loi exponentielle. On le simule exactement ainsi : on additionne des attentes exponentielles tant qu'on reste dans l'intervalle, et on compte les arrivées.

import random
from math import exp, factorial

def nombre_arrivees(lam, duree=1.0):
    """Nombre d'arrivées dans [0 ; duree] quand les attentes suivent la loi exponentielle de paramètre lam."""
    t = random.expovariate(lam)
    compte = 0
    while t <= duree:
        compte += 1
        t += random.expovariate(lam)
    return compte

def proba_poisson(lam, k):
    return exp(-lam) * lam**k / factorial(k)

N = 100_000
minutes = [nombre_arrivees(3) for _ in range(N)]        # 3 requêtes par minute en moyenne
print(sum(minutes) / N)                                 # proche de 3
for k in range(8):
    print(k, round(minutes.count(k) / N, 4), round(proba_poisson(3, k), 4))

Les fréquences simulées collent aux probabilités de la loi de Poisson P(3) : P(X = 0) ≈ 0,0498, P(X = 3) ≈ 0,2240. C'est la vérification par l'expérience du lien entre les deux lois, expliqué au chapitre 8.

Simuler une loi normale avec des uniformes

Le programme propose d'introduire la loi normale en observant la somme de plusieurs nombres de loi uniforme. Une somme de 12 nombres uniformes sur [0 ; 1] a pour espérance 12 × 0,5 = 6 et pour variance 12 × 1/12 = 1 ; si on retire 6, on obtient une variable d'espérance 0 et d'écart type 1, dont la loi est très proche de N(0 ; 1).

import random

def normale_approchee(mu, sigma):
    z = sum(random.random() for _ in range(12)) - 6   # proche de N(0 ; 1)
    return mu + sigma * z                             # proche de N(mu ; sigma)

N = 100_000
volumes = [normale_approchee(420, 25) for _ in range(N)]
for k in (1, 2, 3):
    dedans = sum(1 for v in volumes if abs(v - 420) <= k * 25) / N
    print(f"à moins de {k} écart(s) type : {dedans:.3f}")   # environ 0.679, 0.955, 0.998

Le programme retrouve presque les trois valeurs du chapitre 9 : environ 0,679, 0,955 et 0,998 au lieu de 0,683, 0,954 et 0,997. Ce petit écart ne vient pas du hasard : la somme de 12 uniformes n'est qu'une approximation de la loi normale (elle ne dépasse jamais μ ± 6σ). Pour une vraie loi normale, Python fournit random.gauss(mu, sigma) ; la somme de 12 uniformes montre d'où vient la courbe en cloche.

Simuler une probabilité conditionnelle

Une simulation retrouve aussi le résultat surprenant du chapitre 3 : seulement 59,5 % des alertes signalent une vraie attaque. On simule des tentatives, on ne garde que celles qui déclenchent une alerte, et on compte la part de tentatives malveillantes parmi elles : c'est exactement le sens de « sachant ».

import random

random.seed(1)
alertes = 0
alertes_justifiees = 0
for _ in range(1_000_000):
    malveillante = random.random() < 0.03
    if malveillante:
        alerte = random.random() < 0.95
    else:
        alerte = random.random() < 0.02
    if alerte:
        alertes += 1
        if malveillante:
            alertes_justifiees += 1

print(alertes / 1_000_000)          # proche de P(A) = 0,0479
print(alertes_justifiees / alertes)  # proche de P_A(M) ≈ 0,595

Les calculs exacts avec la bibliothèque standard

Pour des résultats exacts, quelques fonctions suffisent. Range-les dans un fichier lois.py et importe-les quand tu en as besoin :

from math import comb, exp, factorial
from statistics import NormalDist

def binom_fdp(n, p, k):
    """P(X = k), X de loi B(n ; p)."""
    return comb(n, k) * p**k * (1 - p) ** (n - k)

def binom_frep(n, p, k):
    """P(X <= k), X de loi B(n ; p)."""
    return sum(binom_fdp(n, p, i) for i in range(k + 1))

def poisson_fdp(lam, k):
    """P(X = k), X de loi de Poisson de paramètre lam."""
    return exp(-lam) * lam**k / factorial(k)

def poisson_frep(lam, k):
    """P(X <= k), X de loi de Poisson de paramètre lam."""
    return sum(poisson_fdp(lam, i) for i in range(k + 1))

def expo_frep(lam, t):
    """P(T <= t), T de loi exponentielle de paramètre lam."""
    return 1 - exp(-lam * t)

print(round(binom_frep(30, 0.2, 4), 4))                  # 0.2552
print(round(1 - poisson_frep(3, 7), 4))                  # 0.0119
print(round(expo_frep(1e-5, 8760), 4))                   # 0.0839
print(round(NormalDist(420, 25).inv_cdf(0.99), 1))       # 478.2

Les noms reprennent ceux de la calculatrice TI (binomFdp, binomFRép) pour que tu passes facilement de l'une à l'autre. La bibliothèque SciPy fournit les mêmes fonctions (scipy.stats.binom.cdf…), mais elle s'installe à part : pour les exercices du BTS, la bibliothèque standard suffit.

Au tableur

LibreOffice Calc et Excel ont les mêmes fonctions, sous les mêmes noms français. Le dernier argument vaut 0 pour une probabilité ponctuelle (ou une densité) et 1 pour une probabilité cumulée.

Les versions récentes d'Excel proposent aussi des noms avec le suffixe .N (LOI.BINOMIALE.N) ; les noms ci-dessus restent acceptés par les deux logiciels.

Pour représenter une loi, comme le demande le programme, saisis les valeurs k de 0 à 15 en colonne A (de A2 à A17), =LOI.BINOMIALE(A2;30;0,2;0) en B2, recopie vers le bas, puis insère un diagramme en colonnes (Insertion > Diagramme, type Colonne, dans LibreOffice ; Insertion > Graphiques > Histogramme groupé dans Excel) : tu retrouves le diagramme en bâtons du chapitre 5. Avec =LOI.POISSON(A2;3;0), tu obtiens celui de la loi de Poisson de paramètre 3.

Pour simuler au tableur, recopie une formule sur des milliers de lignes : =SI(ALEA()<0,05;1;0) simule un paquet perdu avec la probabilité 0,05, et =-LN(1-ALEA())/0,00001 une durée de vie exponentielle. La touche F9 relance tous les tirages ; une moyenne en bas de colonne montre la fréquence ou la moyenne simulée.

Exercices

import random

def au_moins_un(n, p, essais):
    compte = 0
    for _ in range(essais):
        if any(random.random() < p for _ in range(n)):
            compte += 1
    return compte / essais

print(au_moins_un(30, 0.2, 100_000))   # environ 0.9988
print(1 - 0.8**30)                     # valeur exacte : 0.99876...

import random
from math import exp

durees = [random.expovariate(0.5) for _ in range(200_000)]
survivantes = [d for d in durees if d > 3]
print(sum(1 for d in survivantes if d > 5) / len(survivantes))   # environ 0.368
print(exp(-1))                                                   # 0.3679

import random

valeurs = [0, 1, 2, 3, 4]
probas = [0.30, 0.35, 0.20, 0.10, 0.05]
N = 200_000
chargees = sum(1 for _ in range(N)
               if sum(random.choices(valeurs, weights=probas, k=8)) >= 15)
print(chargees / N)   # environ 0.086