Apprendre à réaliser des simulations aléatoires en Python

La fonction random.seed(a=None, version=2) permet d'initialiser le générateur de nombres aléatoires avec une valeur appelée la graine (seed en anglais) :

>>> import random
>>> 
>>> # initialise le générateur en utilisant l'heure système ou un générateur aléatoire physique
>>> random.seed()

Si a est omis ou égal à None, l'heure système actuelle est utilisée. Si des sources aléatoires sont fournies par le système d'exploitation, elles sont utilisées à la place de l'heure système (voir la fonction os.urandom() pour les détails sur la disponibilité). Si a est un entier, il est utilisé directement.

La graine est souvent générée à partir de l'état du système. Par exemple, en lisant la valeur de l'horloge, ou bien d'un générateur de nombres aléatoires physique (les mouvements de la souris, ou un bruit électronique).

Presque toutes les fonctions du module dépendent de la fonction de base random(), qui génère un nombre aléatoire à virgule flottante entre 0 et 1 inclus.

Elle se base sur l'algorithme de Mersenne Twister.

La fonction peut être appelée ainsi :

>>> # tire au hasard un nombre à virgule flottante entre 0 et 1 inclus.
>>> random.random()
0.5317718825456861

Comme on le verra à la fin, la fonction génère, à partir d'une valeur particulière, toujours la même suite de nombres. Par conséquent, les résultats peuvent être prévus par toute personne connaissant la valeur de départ et l'algorithme utilisé pour créer ces nombres.

L'appel à randrange(start, stop[, step]) génère au hasard un nombre entier choisi dans range(start, stop, step):

>>> # tire un entier au hasard dans la séquence 0, 2, 4, 6, 8
>>> random.randrange(0, 10, 2)
4

Cela équivaut à peu près à choice(range(start, stop, step))mais accepte des séquences d'entiers étendues et est optimisée pour les cas courants.

L'appel à randint(a, b) génère au hasard un nombre entier entre a et b inclus :

>>> # tire au hasard un entier compris entre 1 et 10 inclus
>>> random.randint(1, 10)
7

random.choice(seq) renvoie un élément aléatoire de la séquence non vide seq. Si seq est vide, lève IndexError :

>>> seq = ['a','b','c','d','e','f']
>>> random.choice(seq)
'd'

random.choices(population, weights=None, *, cum_weights=None, k=1) renvoie une liste de k éléments choisis avec remise dans population. Si population est vide, lève IndexError :

>>> seq = ['a','b','c','d','e','f']
>>> random.choices(population=seq, k=10) # tirage avec remise de 10 éléments dans seq
['e', 'b', 'b', 'f', 'e', 'a', 'a', 'b', 'f', 'c']

Mélange la séquence x sans créer de nouvelle instance (« sur place ») :

>>> x = ['a','b','c','d','e','f']
>>> random.shuffle(x)
['c', 'e', 'f', 'b', 'a', 'd']

On réalise ici une permutation au hasard des éléments de x.

Notez que même pour des valeurs de len(x) relativement petites, le nombre total de permutations de x peut rapidement devenir plus grand que la période de la plupart des générateurs de nombres aléatoires. Cela implique que la plupart des permutations d'une longue séquence ne peuvent jamais être générées. Par exemple, une séquence de longueur 2080 est la plus grande qui puisse tenir dans la période du générateur de nombres aléatoires Mersenne Twister.

random.sample(population, k, *, counts=None) renvoie une liste de k éléments choisis dans population sans remise. Si population est vide, lève IndexError.

La fonction peut être appelée comme ceci :

>>> seq = ['a','b','c','d','e','f']
>>> random.sample(population=seq, k=5) # tirage sans remise de 5 éléments dans seq
['e', 'c', 'a', 'd', 'f']

La fonction random.gauss(mu=0.0, sigma=1.0) permet de simuler une variable de Gauss d'espérance mu et d'écart-type sigma.

Le code suivant génére une séquence de 100 valeurs suivant une loi normale :

import random
			
# moyenne μ et écart-type σ de la loi normale
m = 40; s = 4.5
	
# génère une séquence de 100 valeurs suivant une loi normale N(m,s^2)
x = [random.gauss(mu=m, sigma=s) for _ in range(100)]

La fonction choice() va permettre d'effectuer un tirage aléatoire dans la liste ['Pile', 'Face'] simulant ainsi un jeu de pile ou face :

# simulation d'un tirage de pile ou face
resultat = random.choice(["Pile","Face"])

Dans notre exemple, on effectue un tirage aléatoire de pile ou face tant que la variable tirage est égale à 'o' et en affichant à chaque fois les scores des joueurs :

import random

print("Simuler un jeu de pile ou face..\n")

# initialisation des compteurs : compteur joueur n°1, compteur joueur n°2, compteur de tirages
cpt_joueur1 = 0; cpt_joueur2 = 0; cpt_tirages = 0

# saisie des noms des joueurs
nom_joueur1 = input("Nom du joueur n°1 (Pile) :")
nom_joueur2 = input("Nom du joueur n°2 (Face) :")

print()

# tirage ('o': oui, 'n': non)
tirage = input("Commencer le jeu (o/n) :")

print()

# initialisation du générateur de nombres pseudo-aléatoires
random.seed()
 
# réalisation de tirages successifs de pile ou face tant que tirage = "o"
while tirage=="o":
    
    # incrémentation du compteur de tirages
    cpt_tirages += 1 

    print("-- Tirage n°{0} --\n".format(cpt_tirages))
    
    # simulation du tirage de pile ou face
    resultat = random.choice(["Pile","Face"])

    print("Résultat :", resultat)

    # si pile
    if resultat=="Pile":
        cpt_joueur1 += 1 # incrémente le compteur du joueur n°1
    else: # sinon    
        cpt_joueur2 += 1 # incrémente le compteur du joueur n°2

    print()

    # affiche les scores    
    print("{0} : {1} points".format(nom_joueur1, cpt_joueur1))
    print("{0} : {1} points".format(nom_joueur2, cpt_joueur2))

    print()

    # poursuit-on le jeu ('o': oui, 'n': non)
    tirage = input("Poursuivre le jeu (o/n) :")
    
    print()

# affiche le résultat du jeu
if cpt_joueur1>cpt_joueur2: # joueur n°1 gagne
    print('{0} gagne avec {1} points sur {2} possibles !'.format(nom_joueur1, cpt_joueur1, cpt_tirages))
else: # sinon
    if cpt_joueur1<cpt_joueur2: # joueur n°2 gagne
        print('{0} gagne avec {1} points sur {2} possibles !'.format(nom_joueur2, cpt_joueur2, cpt_tirages))
    else: # égalité
        print('Egalité entre les deux joueurs : {0} points chacun !'.format(cpt_joueur1))

L'exécution du code affiche :

Simuler un jeu de pile ou face..

Nom du joueur n°1 (Pile) :Florent
Nom du joueur n°2 (Face) :Olivier

Commencer le jeu (o/n) :o

-- Tirage n°1 --

Résultat : Pile

Florent : 1 points
Olivier : 0 points

Poursuivre le jeu (o/n) :o

-- Tirage n°2 --

Résultat : Pile

Florent : 2 points
Olivier : 0 points

Poursuivre le jeu (o/n) :o

-- Tirage n°3 --

Résultat : Face

Florent : 2 points
Olivier : 1 points

Poursuivre le jeu (o/n) :o

-- Tirage n°4 --

Résultat : Pile

Florent : 3 points
Olivier : 1 points

Poursuivre le jeu (o/n) :n

Florent gagne avec 3 points sur 4 possibles !

La fonction sample() va donc permettre d'effectuer un tirage sans remise de n valeurs dans population :

# tirage aléatoire d'un échantillon de taille n dans population
echantillon = random.sample(population, k=n)

On va maintenant simuler le tirage d'un échantillon dans la population des français afin d'estimer leur taille moyenne :

import random
import numpy as np
import math

print("Simuler un échantillonnage aléatoire..\n")
					
# taille de l'échantillon
n = 1000
 
print("Taille de l'échantillon, n =", n)

# taille moyenne des français et écart-type (valeurs inconnues utilisées pour la simulation)
μ = 175
σ = 7.2

# taille de la population
k = 45000000

print()

print("Création de la population de base de taille k = {0} ..\n".format(k)) 
# création de la  population de base : tailles des français
# en supposant que leur taille se répartisse suivant une loi normale de moyenne m et d'écart type s.
# population = [random.gauss(mu=m, sigma=s) for _ in range(k)]
population = list(np.random.normal(loc=μ, scale=σ, size=k))
 
# tirage aléatoire d'un échantillon de dimension n dans population
echantillon = random.sample(population, k=n)
 
print("Echantillon représentatif :")
 
# Echantillon représentatif
print(echantillon)

print()

# calcul de la moyenne et de l'écart-type de l'échantillon
m = np.mean(echantillon)
s = np.std(echantillon)

print("Moyenne de l'échantillon =", m)

# intervalle de confiance sur la moyenne à 99,7 %
print("Intervalle de confiance sur la moyenne =", [m - 3*σ/math.sqrt(n), m + 3*σ/math.sqrt(n)])

print()

print("Ecart-type de l'échantillon =", s)

En supposant que la taille des français se répartisse suivant une loi normale, on utilise ici la fonction random.normal(loc=0.0, scale=1.0, size=None) du module numpy pour générer la population de base.

L'intervalle kitxmlcodeinlinelatexdvp[m-3\frac{\sigma}{\sqrt{n}};m+3\frac{\sigma}{\sqrt{n}}]finkitxmlcodeinlinelatexdvp est un intervalle de confiance de la moyenne m à environ 99,7 %. Cela veut dire qu'il y a environ 99,7 % de chances que la vraie valeur de la moyenne soit dans cet intervalle.

L'exécution du code affiche :

Simuler un échantillonnage aléatoire..
					
Taille de l'échantillon, n = 1000

Création de la population de base de taille k = 45000000 ..

Echantillon représentatif :
[177.7548635176891, 172.95797319858872, 189.45775722849214, 173.29714340318452, 182.44108042734678, 180.090200994285, 187.1911861793668, 182.2286639672886, 164.66804791571886, 187.23141106882252,...]

Moyenne de l'échantillon = 175.00715275211775
Intervalle de confiance sur la moyenne = [174.32410077752138, 175.69020472671411]

Ecart-type de l'échantillon = 7.24380196890074

Le principe est d'utiliser la fonction gauss() du module random afin d'obtenir une séquence de valeurs suivant une loi normale :

import random
...			
# nombre de tirages successifs dans la simulation
k = 100000

# tirage de k valeurs suivant une loi normale
for i in range(k):
    # génère au hasard une valeur suivant la loi normale N(μ,σ^2)
    xi = random.gauss(μ,σ)
	...

L'objectif va être de générer un graphique en Python permettant de montrer que les valeurs obtenues à l'aide de la fonction gauss() se répartissent bien suivant une courbe de Gauss.

Déroulé du code mise en œuvre :

1. Initialisation des paramètres μ et σ ;
2. Création de la liste des valeurs en abscisse et initialisation des valeurs en ordonnées ;
3. Tirage de k nombres réels suivant une loi normale N(μ,σ²) (k très grand) ;
4. Traçage du graphique représentant la distribution de Gauss et les fréquences d'apparition de chaque valeur de x.

ce qui donne :

import random
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import norm

print("Simuler une gaussienne..\n")
					
# nombre d'épreuves de Bernoulli
n = 80

# probabilité de succès
p = 0.5

# moyenne μ et écart-type σ de la loi normale : N(μ,σ^2)≃ B(n,p)
μ = n*p
σ = math.sqrt(n*p*(1-p))

print("Espérance, μ =", μ)
print("Écart-type, σ =", σ)

# écart entre deux valeurs consécutives
h = 0.5

# génère la séquence des x : 0 -> n avec un pas de h
x = np.arange(0, n+1, h)
y = np.array([0]*len(x)) # [0, ..., 0]

print()
print("Tirage des valeurs suivant la loi normale..")

# nombre de tirages successifs : simulation d'une variable de Gauss
k = 1000000

# initialisation du générateur de nombres pseudo-aléatoires
random.seed()

# tirage de k valeurs suivant une loi normale N(μ,σ^2)
for _ in range(k):
    # génère au hasard une valeur suivant la loi normale N(μ,σ^2)
    xi = random.gauss(μ,σ)
    
    # détermination de l'indice de la valeur dans la liste x
    i = round(xi/h) # x[i] - h/2 ≤ xi ≤ x[i] + h/2

    # incrémentation du compteur de la valeur d'indice i dans y
    y[i] += 1

# calcul de la moyenne des valeurs en tenant compte de leur poids
m = np.average(x, weights=y)
 
# calcul de la variance des valeurs en tenant compte de leur poids
v = np.average((x-m)**2, weights=y)
 
# calcul de l'écart-type
s = math.sqrt(v)
 
print()
print("Moyenne observée =", m)
print("Écart-type observé =", s)

print()
print("Génération du graphique..")

# division des effectifs par le nombre k de tirages
y = y/k

# nom et dimension de la figure contenant le graphique
plt.figure(num="Figure : Loi normale", figsize=(8, 5), dpi=80)

# trace la ccurbe de Gauss
plt.plot(x, h*norm.pdf(x,μ,σ), color="red", label='Loi normale N(μ,σ' + chr(0x00B2) + ')')

# trace les barres correspondant aux fréquences observées
plt.bar(x, y, width = h, color="white", edgecolor = 'blue', label='fréquences observées')

# limites sur l'axe des x : [μ - 4σ, μ + 4σ]
plt.xlim(μ - 4*σ, μ + 4*σ)

# limite inférieure sur l'axe des y
plt.ylim(0)

plt.title("Simulation d'une variable de Gauss")

# ajout de la légende
plt.legend()

# affiche le graphique
plt.show()

L'exécution du code affiche :

Simuler une gaussienne..

Espérance, μ = 40.0
Écart-type, σ = 4.47213595499958

Tirage des valeurs suivant la loi normale..

Moyenne observée = 40.000915
Écart-type observé = 4.478495357011661

Génération du graphique..

Puis le graphique :

Ces algorithmes génèrent, à partir d'une ou de plusieurs valeurs initiales et souvent à l'aide d'une relation de récurrence, toujours la même suite de nombres. Chaque terme produit peut être associé à un état du générateur (valeur du terme et éventuellement de termes précédents, etc.) qui dépend entièrement du précédent. Si un état apparaît une deuxième fois, toute la suite se reproduit à partir de cet état. La période est alors égale au nombre de valeurs générées entre entre deux états identiques.

Il est évident que l'on va chercher à obtenir une suite de nombres qui se répartissent toujours bien suivant une loi uniforme et ceci sur une période la plus longue possible. Par exemple, l'algorithme de Mersenne Twister réputé pour sa qualité possède une période constante égale à 2¹⁹⁹³⁷-1.

La qualité des sorties augmente aussi souvent au détriment de la vitesse de génération et ces paramètres doivent également être considérés lors de l'utilisation d'un générateur.

Dans notre cas, on va choisir de créer un générateur congruentiel linéaire car il est facile à mettre en œuvre et généralement de bonne qualité sous réserve bien sûr de choisir les bons paramètres.

D'arès Wikipedia, un générateur congruentiel linéaire (LCG : linear congruential generator en anglais) est un générateur de nombres pseudo-aléatoires dont l'algorithme, introduit en 1948 par Derrick Lehmer, est basé sur des congruences et une fonction affine.

Les nombres pseudo aléatoires forment une suite dont chaque terme dépend du précédent, selon la formule : kitxmlcodeinlinelatexdvpX_{n+1}=(aX_{n}+c)\mod mfinkitxmlcodeinlinelatexdvp

où a est le multiplicateur, c l'incrément et m le module.

Le terme initial s est appelé la graine (seed en anglais). C'est elle qui va permettre de générer une suite apparemment aléatoire. Pour chaque graine, on aura une nouvelle suite. Cependant, il est possible que certaines graines permettent d'obtenir une suite plus aléatoire que d'autres.

Les termes de cette suite sont donc compris entre 0 et m−1. De plus, comme chaque terme dépend entièrement du précédent, si un nombre apparaît une deuxième fois, toute la suite se reproduit à partir de ce nombre. La période, c'est à dire le nombre de valeurs que produit le générateur avant de recommencer la même séquence, est donc dans ce cas au maximum de m.

Le module est souvent choisi de manière à être de l’ordre de la longueur des mots sur l’ordinateur (par exemple : 2³² sur une machine 32 bits). En effet, les ordinateurs calculant naturellement en base binaire, cela évite une division euclidienne dans le calcul du modulo (opération souvent coûteuse en temps de calcul).

La période est donc relativement courte mais cette méthode présente un avantage : on connaît les critères sur les nombres a, c et m qui vont permettre d’obtenir une période maximale (égale à m).

Pour plus de détail vous pouvez consulter cette page WikipediaGénérateur congruentiel linéaire.

Notre générateur congruentiel linéaire est défini à l'aide d'une classe Random  :

III-D-3-a. Méthode constructeur __init__()▲

class Random():
    # classe permettant de définir le générateur congruentiel linéaire (LCG)
    
    def __init__(self, m=2147483648, a=1103515245, c=12345, s=int(time.time())):
        # m: module ; a: multiplicateur; c: incrément; s: valeur initiale ou graine (seed)
        # Si les arguments de la fonction sont omis, on passe les mêmes valeurs que pour le générateur d'UNIX

        # mise à jour des attributs avec les valeurs passées en argument
        self.m = m
        self.a = a
        self.c = c
        self.x = s
)

III-D-3-b. Méthode rand()▲

class Random():
    # classe permettant de définir le générateur congruentiel linéaire (LCG)
    ....
    
    def rand(self):
        # méthode permettant de générer un nombre aléatoire entre 0 et 1
        
        # formule de récurrence pour le LCG : Xn+1 = (aXn + c) mod m
        self.x = (self.a * self.x + self.c) % self.m

        # renvoi la nouvelle valeur pour x divisée par le module m
        return self.x/self.m

Une méthode de Monte-Carlo consiste à calculer une valeur numérique approchée en utilisant un procédé aléatoire. Elle est particulièrement utilisée pour le calcul de surfaces ou de volumes.

Soit un cercle de centre 0 et de rayon 1, inscrit dans un carré de côté 2.

On considère maintenant un quart de ce cercle inscrit dans un petit carré de côté 1.

Le procédé aléatoire va consister à tirer au hasard (et successivement) un grand nombre de points dans le petit carré et à compter le nombre de points situés dans le quart de cercle :

En faisant le rapport du nombre de points situés dans le cercle sur le nombre total de tirages, on doit obtenir une valeur approchée du rapport entre la surface du cercle et celle du carré, c'est à dire Pi/4.

Pour avoir plus de détail vous pouvez consulter cette page WikipediaMéthode de Monte-Carlo.

III-D-4-a. Fonction Test▲

def test_rand(n):
    # n : nombre de tirages

    # création de l'objet Random avec les paramètres par défaut : initialisation du générateur de nombres pseudo-aléatoires
    rand_lcg = Random()

    # initialisation du compteur et de la liste des points
    cpt = 0; points=[]
        
    # réalisation d'un tirage de n points dans le carré de côté 1
    for i in range(n):

        # tirage au hasard de la coordonnée x
        x = rand_lcg.rand()
        
        # tirage au hasard de la coordonnée y
        y = rand_lcg.rand()

        # ajout du point à la liste
        points.append((x,y))
        
        # test si le point de coordonnées (x,y) est dans le cercle
        if x*x + y*y <= 1:
           # si oui on incrémente le compteur 
           cpt +=  1 

    # renvoie la valeur approximation de Pi : (cpt / n) * 4, et la liste des points tirés au hasard
    return ((cpt / n) * 4, points)

III-D-4-b. Simulation de Monte-Carlo▲

On va maintenant tester notre générateur à l'aide de la méthode de Monte-Carlo décrite précédemment et aussi visualiser sur un graphique la répartition uniforme des points sur le carré :

 

Sélectionnez

import math
import numpy as np
import matplotlib.pyplot as plt
						
# 3 tests successifs de notre générateur en faisant varier n
for i in range(3, 6):

    # valeur de n : 10^(i)
    n = pow(10,i)

    print()

    print("Test n°{0} :".format(i-2)) 
    print("Nombre de points choisis au hasard =", n)

    # calcul approché de Pi à l'aide d'une méthode de Monte-Carlo.
    pi, points = test_rand(n)
 
    print("Approximation de Pi =", pi)
    print("Erreur absolue =", abs(pi-math.pi))

    # tri des points en fonction des x
    points.sort()

    # créé la liste des x et des y à partir des points
    x = [p[0] for p in points]
    y = [p[1] for p in points]

    # calcul des y pour tracer le quart de cercle
    fx = [math.sqrt(1-xi*xi) for xi in x]

    # nom et dimension de la figure contenant le graphique
    plt.figure(num="Simulation de Monte-Carlo : calcul de Pi", figsize=(7, 6), dpi=80) 
    
    # même échelle sur l'axe des x et des y
    plt.axis('scaled')

    # trace le quart de cercle en rouge
    plt.plot(x, fx, color="red")

    # trace les points tirés au hasard
    plt.scatter(x, y, s=0.5)

    # limites inférieures et supérieures sur l'axe des x et des y
    plt.xlim(0, 1); plt.ylim(0, 1)

    # affiche le graphique
    plt.show()

Le code affiche :

 

Sélectionnez

Test n°1 :
Nombre de points choisis au hasard = 1000
Approximation de Pi = 3.196
Erreur absolue = 0.05440734641020706

Puis :

 

Sélectionnez

Test n°2 :
Nombre de points choisis au hasard = 10000
Approximation de Pi = 3.162
Erreur absolue = 0.020407346410206806

Et enfin :

 

Sélectionnez

Test n°3 :
Nombre de points choisis au hasard = 100000
Approximation de Pi = 3.14572
Erreur absolue = 0.004127346410206734

A l'exécution on remarque que la fonction tend bien vers Pi à mesure que le nombre de tirages augmente. On peut aussi constater sur les différents graphiques que les points se répartissent à peu près uniformément sur le carré de côté 1. J'ai pensé qu'il serait intéressant d'ajouter à l'estimation de Monte-carlo (calcul approché de Pi) une représentation graphique des points sur le carré afin de mieux visualiser leur répartition uniforme. Cependant, je suis aussi conscient des limites de cette méthode. En réalité, il faudrait procéder à une analyse statistiques prenant en compte d'autres paramètres, comme sa période. Ceci permettrait de déceler certains défauts du générateur.