Simuler une loi de Poisson

import numpy.random as rd
import matplotlib.pyplot as plt
from math import exp, factorial

def poisson(lam):
    # methode de Knuth : on multiplie des uniformes jusqu'a passer sous e^-lam
    seuil, produit, k = exp(-lam), rd.random(), 0
    while produit > seuil:
        produit = produit * rd.random()
        k = k + 1
    return k

N, lam = 100000, 4
effectifs = [0] * 20
for _ in range(N):
    k = poisson(lam)
    if k < 20:
        effectifs[k] += 1

freq = [e / N for e in effectifs]
theo = [exp(-lam) * lam**k / factorial(k) for k in range(20)]
plt.bar([k - 0.2 for k in range(20)], freq, width=0.4, label="simulees")
plt.bar([k + 0.2 for k in range(20)], theo, width=0.4, label="theoriques")
plt.xlabel("k"); plt.legend(); plt.show()

Ce qu'on voit : deux profils qui se superposent, avec deux barres maximales à égalité en $k=3$ et $k=4$ (toutes deux à $0,195$) — car pour $\lambda$ entier, $P(X=\lambda -1)=P(X=\lambda )$ exactement.

L'écart attendu : $\sqrt{p(1-p)/N}\approx 1,3\times 10^{-3}$ pour $p\approx 0,195$. Les différences constatées sont de cet ordre.

Pourquoi la méthode de Knuth fonctionne : le produit de variables uniformes descend sous ${\mathrm{e}}^{-\lambda }$ après un nombre d'étapes qui suit précisément une loi de Poisson. On ne demande jamais à l'élève de la démontrer, mais il faut savoir qu'aucune simulation « naïve » n'est possible : le support est infini.