Simuler une loi géométrique

Le code.

import numpy.random as rd

def geometrique(p):
    # rang du premier succes : on compte jusqu'a reussir
    k = 1
    while rd.random() >= p:     # echec : on recommence
        k = k + 1
    return k

N = 100000
p = 0.2
X = [geometrique(p) for _ in range(N)]
print(sum(X) / N, 1 / p)                  # ~5.0   5.0

Y = rd.geometric(p, N)                    # la meme loi, en une ligne
print(sum(Y) / N)                         # ~5.0

Pourquoi ce code simule bien $\mathcal{G}(p)$. L'appel rd.random() rend un réel tiré uniformément dans $[0,1[$ ; l'événement « ce réel est $<p$ » a donc la probabilité $p$. Chaque tour de boucle est ainsi une épreuve de Bernoulli de paramètre $p$, indépendante des précédentes, et la variable k compte les épreuves jusqu'au premier succès inclus : c'est la définition même du rang d'apparition du premier succès.

Ce qu'on attend. Le cours donne $E(X)=\frac{1}{p}=5$ et $V(X)=\frac{1-p}{p^2}=\frac{0,8}{0,04}=20$, donc $\sigma (X)=\sqrt{20}\approx 4,47$. La moyenne empirique de $N=10^5$ tirages a un écart-type de l'ordre de $\frac{\sigma (X)}{\sqrt{N}}=\frac{4,47}{316}\approx 0,014$ : on doit donc lire une valeur comprise entre $4,96$ et $5,04$ environ. Un écart de cet ordre ne signale aucune erreur.

Le test $\ge p$ plutôt que $>p$. L'événement « rd.random() vaut exactement $p$ » a une probabilité nulle en théorie : les deux écritures donnent la même loi. Ce qui compte est que l'échec soit le complémentaire exact du succès, sans quoi la boucle ne compterait plus la bonne chose.

La boucle se termine-t-elle ? À chaque tour, la probabilité de sortir vaut $p>0$ ; la probabilité de faire encore $n$ tours est $(1-p{)}^n$, qui tend vers $0$. L'arrêt n'est donc pas garanti au sens strict, mais il est presque sûr — la nuance à laquelle une simulation de loi géométrique oblige.