Échantillonnage et simulation

Comment un institut de sondage annonce-t-il le salaire moyen des Français en n'interrogeant que quelques milliers de personnes ? Comment un fabricant contrôle-t-il le poids moyen de ses paquets de café sans peser toute la production ? Comment vérifier qu'un dé est équilibré sans démonstration possible, seulement en le lançant ? Ces trois questions reposent sur la même idée : observer un échantillon pour estimer une grandeur inaccessible. Ce chapitre met cette idée à l'épreuve de l'expérience, grâce à un outil que le hasard ne peut pas fatiguer : la simulation en Python. Au passage, il introduit la seule nouveauté algorithmique de l'année, la structure de liste, indispensable pour stocker et exploiter les résultats de milliers d'expériences aléatoires.

Simuler le hasard avec Python

Lancer un dé $10000$ fois à la main est fastidieux ; demander à un ordinateur de le faire prend quelques millisecondes. Le module random de Python fournit les briques de base de toute simulation.

from random import random

def tir_reussi():
    # une basketteuse reussit ses lancers francs avec probabilite 0.8
    return random() < 0.8

Pour simuler une variable aléatoire quelconque, on découpe l'intervalle $[0;1[$ en zones dont les longueurs sont les probabilités de la loi : c'est la méthode du cumul des probabilités.

from random import random

def gain():
    r = random()
    if r < 0.7:
        return 0
    elif r < 0.95:
        return 2
    else:
        return 10

Lorsque les probabilités sont des fractions simples de même dénominateur, choice offre un raccourci : il suffit de faire apparaître chaque valeur dans la liste proportionnellement à sa probabilité.

from random import choice

def simule_Y():
    return choice([0, 1, 1, 5])

Les listes en Python

Simuler mille lancers de dé produit mille résultats : il faut une structure pour les ranger tous dans une seule variable. C'est le rôle de la liste, suite finie et ordonnée de valeurs, notée entre crochets.

Une liste de longueur $n$ a donc des indices allant de $0$ à $n-1$ : demander L[3] dans la liste précédente provoque une erreur. Deux fonctions accompagnent constamment les listes : len(L) renvoie le nombre d'éléments de L, et sum(L) renvoie la somme de ses éléments (lorsqu'ils sont des nombres).

L = [2, 7, 5]
print(L[0])     # affiche 2
print(len(L))   # affiche 3
print(sum(L))   # affiche 14
L[1] = 10       # modifie le deuxieme element
print(L)        # affiche [2, 10, 5]

Les trois méthodes ci-dessous construisent la même liste des carrés des entiers de $0$ à $9$ :

# en extension (penible au-dela de quelques elements)
carres = [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

# par ajouts successifs
carres = []
for n in range(10):
    carres.append(n**2)

# en comprehension
carres = [n**2 for n in range(10)]

Une condition peut filtrer les éléments retenus :

# carres des entiers pairs entre 0 et 9
carres_pairs = [n**2 for n in range(10) if n % 2 == 0]
print(carres_pairs)  # affiche [0, 4, 16, 64]

notes = [12, 8, 15, 10, 17]

# iteration directe : compter les notes au-dessus de 10
compteur = 0
for note in notes:
    if note > 10:
        compteur = compteur + 1
print(compteur)  # affiche 3

# parcours par indices : afficher chaque note avec sa position
for i in range(len(notes)):
    print("Note numero", i, ":", notes[i])

def moyenne(L):
    return sum(L) / len(L)

print(moyenne([12, 8, 15, 10, 17]))  # affiche 12.4

Échantillons et moyenne d'échantillon

Les outils sont en place : le hasard se simule, les résultats se rangent dans des listes. Reste à définir l'objet central du chapitre.

Attention à ne pas confondre deux nombres de natures différentes : l'espérance $\mu =E(X)$ est un paramètre théorique, calculé une fois pour toutes à partir de la loi de $X$ ; la moyenne $m$ est une valeur observée, qui dépend de l'échantillon obtenu et change d'un échantillon à l'autre.

from random import randint

def moyenne_echantillon(n):
    # simule n lancers d'un de equilibre et renvoie leur moyenne
    echantillon = [randint(1, 6) for k in range(n)]
    return sum(echantillon) / n

print(moyenne_echantillon(100))  # par exemple : 3.42

Fluctuation d'échantillonnage et estimation de l'espérance

L'exemple du dé le montre : deux échantillons de même taille donnent en général deux moyennes différentes. Ce phénomène porte un nom.

Que devient cette fluctuation quand la taille de l'échantillon grandit ? Simulons un très long échantillon de lancers de dé et traçons la moyenne des $n$ premiers lancers en fonction de $n$ :

Au début, la moyenne oscille fortement : elle ne dépend que de quelques lancers. Mais à mesure que $n$ grandit, les oscillations s'amortissent et la courbe se stabilise au voisinage de l'espérance $\mu =\mathrm{3,5}$. En classe de seconde, le même constat avait été fait pour la fréquence d'une issue dans une expérience à deux issues ; il se prolonge ici à la moyenne d'une variable aléatoire quelconque.

Ce résultat fonde le principe d'estimation : pour connaître une espérance, ou une moyenne dans une population trop vaste pour être observée entièrement, on l'estime par la moyenne calculée sur un grand échantillon.

Reste une question, cruciale dans les trois situations : à partir de quel écart entre $m$ et la valeur attendue doit-on s'étonner ? Autrement dit, quelle est l'ampleur normale de la fluctuation d'échantillonnage ? C'est l'objet de la dernière section.

Mesurer la dispersion : l'écart $2\sigma \mathrm{/}\sqrt{n}$

Pour mesurer la fluctuation, une seule moyenne ne suffit pas : il en faut beaucoup. On simule donc $N$ échantillons de taille $n$, on calcule les $N$ moyennes obtenues, et on les range... dans une liste, naturellement.

from random import randint

def moyenne_echantillon(n):
    echantillon = [randint(1, 6) for k in range(n)]
    return sum(echantillon) / n

# liste des moyennes de N = 1000 echantillons de taille n = 50
moyennes = [moyenne_echantillon(50) for k in range(1000)]

L'histogramme de ces $1000$ moyennes révèle une structure remarquable :

Les moyennes s'accumulent autour de $\mu =\mathrm{3,5}$ et se raréfient à mesure qu'on s'en éloigne. Les traits en pointillés délimitent la zone $[\mu -{\displaystyle \frac{2\sigma }{\sqrt{n}}};\mu +{\displaystyle \frac{2\sigma }{\sqrt{n}}}]$, où $\sigma$ désigne l'écart type de $X$ : la grande majorité des moyennes s'y trouve. Calculons ces bornes pour le dé.

Quelle proportion des moyennes tombe dans cette zone ? Écrivons la fonction qui répond à la question, conformément au programme.

from random import randint

def moyenne_echantillon(n):
    echantillon = [randint(1, 6) for k in range(n)]
    return sum(echantillon) / n

def proportion_ecarts(N, n):
    # proportion des N echantillons de taille n dont la moyenne m
    # verifie |m - mu| <= 2*sigma/racine(n), pour un de equilibre
    mu = 3.5
    sigma = (35 / 12)**0.5
    compteur = 0
    for k in range(N):
        m = moyenne_echantillon(n)
        if abs(m - mu) <= 2 * sigma / n**0.5:
            compteur = compteur + 1
    return compteur / N

print(proportion_ecarts(1000, 50))  # par exemple : 0.957

En exécutant ce programme plusieurs fois, avec diverses valeurs de $N$ et de $n$, on obtient des proportions comme $\mathrm{0,957}$, $\mathrm{0,948}$, $\mathrm{0,951}$... Le constat est frappant, et il se reproduit avec d'autres variables aléatoires que le dé.

Ce résultat quantifie la fluctuation d'échantillonnage : l'écart « normal » entre une moyenne observée et l'espérance est de l'ordre de ${\displaystyle \frac{2\sigma }{\sqrt{n}}}$. Deux enseignements pratiques en découlent.

La précision se paie au carré. Le dénominateur $\sqrt{n}$ grandit lentement : pour diviser l'écart ${\displaystyle \frac{2\sigma }{\sqrt{n}}}$ par $2$, il faut multiplier $\sqrt{n}$ par $2$, donc multiplier la taille de l'échantillon par $4$. Pour le dé, l'écart vaut environ $\mathrm{0,483}$ avec $n=50$, environ $\mathrm{0,242}$ avec $n=200$, et il faut $n=800$ pour descendre vers $\mathrm{0,121}$. C'est pourquoi les sondages, coûteux, se contentent souvent de quelques milliers de personnes : gagner un chiffre de précision multiplierait le coût par $100$.

Un critère pour s'étonner. Reprenons le dé suspect de la section précédente : sur $1000$ lancers, l'écart normal vaut ${\displaystyle \frac{2\sigma }{\sqrt{1000}}}\approx {\displaystyle \frac{2\times \mathrm{1,708}}{\mathrm{31,6}}}\approx \mathrm{0,108}$. Une moyenne observée de $\mathrm{4,2}$, soit un écart de $\mathrm{0,7}$ par rapport à $\mathrm{3,5}$, est plus de six fois supérieure à cet écart : un tel échantillon est exceptionnel pour un dé équilibré. La conclusion raisonnable est que le dé n'est pas équilibré.

def moyenne_gains(n):
    echantillon = [gain() for k in range(n)]
    return sum(echantillon) / n

print(moyenne_gains(10000))  # par exemple : 1.0134

L'essentiel en un coup d'œil

• Simuler : random() (réel uniforme dans $[0;1[$), randint(a, b) (entier entre $a$ et $b$ inclus), choice(liste) ; une expérience de probabilité de succès $p$ se simule par le test random() < p ; une loi quelconque, par cumul des probabilités ou par choice sur une liste aux valeurs répétées proportionnellement à leurs probabilités.
• Listes : génération en extension [1, 2, 3], par append dans une boucle, ou en compréhension [f(k) for k in range(n)] (avec condition éventuelle) ; indices à partir de $0$ ; len, sum ; parcours par indices ou itération directe sur les éléments.
• Échantillon de taille $n$ : liste des $n$ résultats de $n$ répétitions indépendantes de l'expérience, dans les mêmes conditions ; sa moyenne $m$ varie d'un échantillon à l'autre (fluctuation d'échantillonnage).
• Estimation : lorsque $n$ est grand, sauf exception, $m$ est proche de l'espérance $\mu$ — on estime donc $\mu$ (ou une moyenne inconnue dans une population) par la moyenne observée sur un grand échantillon.
• Dispersion : en simulant $N$ échantillons de taille $n$, environ $95\mathrm{\%}$ des moyennes vérifient $\mathrm{\mid }m-\mu \mathrm{\mid }⩽{\displaystyle \frac{2\sigma }{\sqrt{n}}}$ (constat expérimental) ; diviser cet écart par $2$ exige un échantillon $4$ fois plus grand.