Table des matières PDFPython

Échantillonnage des distributions de probabilité discrètes

1. Introduction

On s'intéresse à la génération de nombres pseuso-aléatoires suivant une distribution de probabilité non uniforme. On dispose d'un générateur de nombres pseudo-aléatoires qui délivre une séquence de nombres entiers compris entre 0 et M, avec une distribution uniforme. On utilisera le générateur du module python random. Ces nombres entiers sont convertis (par division par M) en nombres à virgule flottante compris entre 0 et 1 avec la fonction random.random(), ou bien en nombres entiers sur un intervalle quelconque avec la fonction random.randint(a,b).

On considère N évènements référencés par un indice k variant de 0 à N-1. La probabilité de l'évènement k est noté pk et on a :

On cherche à tirer des nombres entiers aléatoires compris entre 0 et N-1 avec la distribution de probabilité (pk).

Nous allons voir trois méthodes pour y parvenir : la méthode du rejet, la méthode d'inversion de la fonction de répartition et la méthode de Metropolis.

On verra aussi comment ces méthodes peuvent s'appliquer aux distributions de probabilité continues discrétisées.

2. Méthode du rejet

Soit pmax la plus grande probabilité de la distribution (pk). La figure suivante montre une représentation graphique d'une distribution à 6 évènements.

figureA.svgFigure pleine page

La méthode du rejet consiste à tirer un nombre aléatoire entier compris entre 0 et N-1 avec une probabilité uniforme. Si k est le nombre tiré, un nombre aléatoire réel x est tiré avec une densité de probabilité uniforme dans l'intervalle [0,pmax]. Si x<pk alors le nombre k est accepté, sinon il est rejeté.

Graphiquement, cela revient à lancer des jetons sur le rectangles [0,N-1]x[0,pmax] et à rejeter les jetons qui tombent en dehors des rectangles coloriés. Sur la figure précédente, 7 tirages sont effectués, et seulement 3 sont acceptés (coloriés en rouge).

On remarque qu'il n'est pas nécessaire que la distribution soit normalisée.

from matplotlib.pyplot import *
import time
import math
import random
               

La fonction random de la classe suivante renvoit un nombre aléatoire avec une distribution fournie en argument dans le constructeur, sous forme d'une liste.

class rejet:
    def __init__(self,proba):
    
    def random(self):
               
class rejet:
    def __init__(self,proba):
        self.proba = proba
        self.n = len(proba)
        self.pmax = max(proba)
    def random(self):
        fin = False
        while not fin:
            k = random.randint(0,self.n-1)
            x = random.random()*self.pmax
            fin = x <= self.proba[k]
        return k
               

Pour tester ce générateur, on effectue un million de tirages que l'on place dans une pile afin d'obtenir l'histogramme. On mesure aussi le temps de calcul.

proba = [1,2,3,1,1,2]
random.seed(143652)
liste = []
N = int(1e6)
t = time.clock()
rand = rejet(proba)
for i in range(N):
    k = rand.random()
    liste.append(k)
t = time.clock()-t
figure()
hist(liste,len(proba))
              
print(t)
--> 4.800947
figAfigA.pdf

Cette méthode est efficace lorsque le taux de rejet est faible, c'est-à-dire lorsque la distribution est presque uniforme.

3. Inversion de la fonction de répartition

La fonction de répartition (ou fonction de distribution cumulative) donne la probablité d'obtenir un évènement dans l'intervalle [0,n-1] :

Voyons tout d'abord une approche graphique. Au lieu de disposer des rectangles de hauteur pk côte à côte, on empile horizontalement des rectangles de largeur pk.

figureB.svgFigure pleine page

Des tirages sont faits dans le rectangle (de largeur FN) : un nombre réel x aléatoire est généré dans l'intervalle [0,FN] avec une densité de probabilité uniforme. Il faut déterminer dans quel rectangle le point correspondant est situé, c'est-à-dire trouver l'indice k vérifiant :

L'indice k est le nombre aléatoire obtenu pour ce tirage.

On voit aussi que l'indice k est obtenu par inversion de la fonction de répartition.

Lorsque la taille de la distribution est faible, on peut rechercher l'indice k en parcourant les éléments de la fonction de répartition dans l'ordre. La classe suivante effectue ce calcul :

class inversion:
    def __init__(self,proba):
    
    def random(self):
                  
class inversion:
    def __init__(self,proba):
        self.proba = proba
        self.n = len(proba)
        self.repartition = [0]
        for k in range(1,self.n+1):
            self.repartition.append(self.repartition[k-1]+self.proba[k-1])
        self.max = self.repartition[self.n]
    def random(self):
        x = random.random()*self.max
        k = 0
        while x>=self.repartition[k+1]: 
            k += 1
        return k
                  
liste = []
t = time.clock()
rand = inversion(proba)
for i in range(N):
    k = rand.random()
    liste.append(k)
t = time.clock()-t
figure()
hist(liste,len(proba))
              
print(t)
--> 1.2371030000000003
figBfigB.pdf

On voit sur cet exemple que la méthode d'inversion est beaucoup plus rapide que la méthode du rejet. Lorsque le nombre d'évènements est grand, on améliore la méthode en faisant une recherche de l'indice k par dichotomie :

class inversion_dichotomie:
    def __init__(self,proba):
    
    def random(self):
               
class inversion_dichotomie:
    def __init__(self,proba):
        self.proba = proba
        self.n = len(proba)
        self.repartition = [0]
        for k in range(1,self.n+1):
            self.repartition.append(self.repartition[k-1]+self.proba[k-1])
        self.max = self.repartition[self.n]
    def random(self):
        x = random.random()*self.max
        kmin = 0
        kmax = self.n
        fin = False
        while not fin:
            k = int((kmin+kmax)/2)
            if self.repartition[k+1]<=x:
                kmin = k
            elif x<self.repartition[k]:
                kmax = k
            else:
                fin = True
        return k
                  

Voyons comment utiliser cette méthode pour échantillonner une distribution continue (discrétisée). On considère comme exemple une distribution en cosinus carré. La distribution est discrétisée sur un intervalle donné (100 valeurs sur 2 périodes).

nper = 2
nechant = 100
dt = math.pi*nper/nechant
proba = []
for i in range(nechant):
    proba.append(math.cos(dt*i)**2)
liste = []
t = time.clock()
rand = inversion_dichotomie(proba)
for i in range(N):
    k = rand.random()
    liste.append(k)
t = time.clock()-t
figure()
hist(liste,len(proba))
                  
print(t)
--> 3.5766349999999996
figCfigC.pdf

4. Méthode de Metropolis

La méthode de Metropolis a été introduite pour effectuer des calculs de physique statistique. Voir Algorithme de Metropolis pour une présentation détaillée de la méthode. On se contente ici de donner un exemple d'utilisation pour obtenir l'échantillonnage d'une distribution non uniforme.

La méthode consiste à générer une suite de nombres aléatoires (une chaîne de Markov). Un nombre k' s'obtient à partir du précédent k de la manière suivante :

Il faut remarquer que si le nouveau nombre k' est refusé alors le nombre suivant est identique au précédent. Cela est tout à fait différent de la méthode de rejet, où un nouveau tirage doit être fait à chaque fois qu'une valeur est rejetée.

La classe suivante effectue ce calcul. La première valeur de la suite est la moitié de N.

class metropolis:
    def __init__(self,proba):
        self.proba = proba
        self.n = len(proba)
        self.indice = int(self.n/2)
    def random(self):
        i = random.randint(0,self.n-1)
        if self.proba[i] >= self.proba[self.indice]:
            self.indice = i
        else:
            x = random.random()
            if x < self.proba[i]/self.proba[self.indice]:
                self.indice = i
        return self.indice
               

Voyons l'application à la distribution en cosinus carré :

liste = []
t = time.clock()
rand = metropolis(proba)
for i in range(N):
    k = rand.random()
    liste.append(k)
t = time.clock()-t
figure()
hist(liste,len(proba))
                  
print(t)
--> 3.224786
figDfigD.pdf

La méthode de Metropolis converge plus lentement que la métode d'inversion. Elle est surtout utilisée lorsque le nombre d'états (le nombre d'évènements) est très grand. Par exemple en physique statistique, le nombre d'états est si grand que la méthode d'inversion est absolument inutilisable.

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.