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

Denis Hulo

Objectif : apprendre à réaliser des simulations aléatoires en Python (simulation de variable aléatoire, méthode de Monte-Carlo, etc.).

Niveau requis :confirmé

Commentez cet article : Commentez

- Titre : Apprendre à réaliser des simulations aléatoires en Python

- Auteur : Denis Hulo

- Parution : 9 mai 2025

- Public visé :

- Licence : Copyright ® 2025 Denis Hulo. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- I. Introduction

- II. Module random

- II-A. Générateur de nombres pseudo-aléatoires

- II-A-1. seed()

- II-A-2. random()

- II-B. Fonctions pour les entiers

- II-B-1. randrange()

- II-B-2. randint()

- II-C. Fonctions pour les séquences

- II-C-1. Fonction choice()

- II-C-2. Fonction choices()

- II-C-3. Fonction shuffle()

- II-C-4. Fonction sample()

- II-D. Distributions pour les nombres réels

- II-D-1. Distribution de Gauss

- III. Mise en œuvre de simulations aléatoires

- III-A. Simuler un jeu de pile ou face

- III-A-1. Utilisation de la fonction choice()

- III-A-2. Simulation du jeu

- III-B. Simuler un échantillonnage aléatoire

- III-B-1. Utilisation de la fonction sample()

- III-B-2. Simulation du tirage

- III-C. Simuler une gaussienne

- III-C-1. Utilisation de la fonction gauss()

- III-C-2. Simulation de la gaussienne

- III-D. Créer et tester son propre générateur de nombres aléatoires

- III-D-1. Les PRNG en général

- III-D-2. Générateur congruentiel linéaire

- III-D-3. Classe Random

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

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

- III-D-4. Test du générateur à l'aide d'une simulation de Monte-Carlo

- III-D-4-a. Fonction Test

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

- IV. Conclusion

- VI. Téléchargement

- VII. Remerciements

I - Introduction

Le module random met en œuvre des générateurs de nombres pseudo-aléatoires pour différentes distributions (loi uniforme, loi normale, etc.).

Comme nous le verrons plus loin, ces suites de nombres sont au moins en partie prévisibles, c'est pourquoi on parle de nombres pseudo-aléatoires.

Après avoir présenté quelques fonctions importantes du module random, notre objectif sera de les utiliser pour réaliser différentes simulations (jeu de pile ou face, échantillonnage aléatoire, etc.).

Pour conclure, nous proposerons de créer et de tester notre propre générateur de nombres pseudo-aléatoires à l'aide d'une méthode de Monte-Carlo.

Par la suite, pour simplifier la lecture, on parlera parfois de nombres aléatoires au lieu de nombres pseudo-aléatoires.

II - Module random

En plus de sa fonction de base random(), ce module contient aussi des méthodes permettant par exemple de tirer au hasard des entiers ou encore de choisir des éléments dans une séquence.

II-A - Générateur de nombres pseudo-aléatoires

On présente maintenant les fonctions de base.

II-A-1 - seed()

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).

II-A-2 - random()

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.

II-B - Fonctions pour les entiers

Ce module contient également des méthodes permettant de générer des entiers au hasard.

II-B-1 - randrange()

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.

II-B-2 - randint()

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

II-C - Fonctions pour les séquences

On présente enfin des fonctions permettant de générer des éléments (nombres, caractères, etc.) choisis au hasard dans des séquences.

II-C-1 - Fonction choice()

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'

II-C-2 - Fonction choices()

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']

II-C-3 - Fonction shuffle()

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.

II-C-4 - Fonction sample()

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']

II-D - Distributions pour les nombres réels

II-D-1 - Distribution de Gauss

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)]

III - Mise en œuvre de simulations aléatoires

III-A - Simuler un jeu de pile ou face

Le pile ou face est un jeu de hasard se jouant avec une pièce de monnaie. Le principe du jeu est de lancer en l'air une pièce équilibrée et de parier sur le côté sorti. La pièce tournoyante tombe au sol et s'y stabilise, ou bien elle est rattrapée d'une main et posée à plat dans l'autre main.

On souhaite utiliser la fonction choice afin de simuler un jeu de pile ou face. Chaque pile rapporte 1 point au joueur n°1 et chaque face 1 point au joueur n°2.

III-A-1 - Utilisation de la fonction choice()

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"])

III-A-2 - Simulation du jeu

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("-- Lancé 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

-- Lancé n°1 --

Résultat : Pile

Florent : 1 points
Olivier : 0 points

Poursuivre le jeu (o/n) :o

-- Lancé n°2 --

Résultat : Pile

Florent : 2 points
Olivier : 0 points

Poursuivre le jeu (o/n) :o

-- Lancé n°3 --

Résultat : Face

Florent : 2 points
Olivier : 1 points

Poursuivre le jeu (o/n) :o

-- Lancé 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 !

III-B - Simuler un échantillonnage aléatoire

On souhaite simuler un tirage sans remise à l'aide de la fonction sample(). Ce type de tirage est très utilisé en statistiques quand il s'agit d'obtenir un échantillon représentatif d'une population plus large.

D'après Wikipedia, en statistique, l'échantillonnage ou le sondage, désigne les méthodes de sélection d'un sous-ensemble d'individus (un échantillon) à l'intérieur d'une population pour estimer les caractéristiques de l'ensemble de la population. Cette méthode présente plusieurs avantages : une étude restreinte sur une partie de la population, un moindre coût, une collecte des données plus rapide que si l'étude avait été réalisé sur l'ensemble de la population, la réalisation de contrôles destructifs, etc.

Les résultats obtenus constituent un échantillon. Sur un échantillon, on peut calculer différents paramètres statistiques de position (moyenne, etc.) ou de dispersion (écart type, etc.) issus de la statistique descriptive, de la même manière que l'on peut déterminer des paramètres statistiques d'une population par son recensement exhaustif.

III-B-1 - Utilisation de la fonction sample()

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)

III-B-2 - Simulation du tirage

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 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

III-C - Simuler une gaussienne

On souhaite maintenant simuler une variable aléatoire suivant une loi normale (ou loi de Laplace-Gauss), puis représenter sur un graphique la répartition des valeurs obtenues.

La loi normale est une loi de probabilité continue qui dépend de deux paramètres : son espérance, un nombre réel noté μ, et son écart-type, un nombre réel positif noté σ. La densité de probabilité de la loi normale d'espérance μ et d'écart-type σ est donnée par :

Parmi les lois de probabilité, les lois normales prennent une place particulière grâce au théorème central limite. En effet, elles correspondent au comportement, sous certaines conditions, d'une suite d'expériences aléatoires similaires et indépendantes lorsque le nombre d'expériences est très élevé. Grâce à cette propriété, une loi normale permet d'approcher d'autres lois comme la loi binomiale et ainsi de modéliser de nombreuses études scientifiques comme des mesures d'erreurs ou des tests statistiques, en utilisant par exemple les tables de la loi normale centrée réduite.

III-C-1 - Utilisation de la fonction gauss()

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(μ,σ)
...

III-C-2 - Simulation de la gaussienne

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(μ,σ^2) (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 :

Le graphique montre bien que les barres correspondant aux fréquences obtenues (barres au contour bleu) s'ajustent très bien avec la courbe représentant la distribution de la loi normale (en rouge). La simulation est donc très bonne pour un grand nombre de tirages

III-D - Créer et tester son propre générateur de nombres aléatoires

On se propose enfin de créer son propre générateur de nombres pseudo-aléatoires, afin notamment de mieux comprendre le principe de fonctionnement d'un PRNG (pseudorandom number generator) avec ces différents paramètres (graine, période, etc.). Ensuite, on va également montrer comment le tester à l'aide d'une simulation de Monte-Carlo.

III-D-1 - Les PRNG en général

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^19937-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.

III-D-2 - Générateur congruentiel linéaire

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 :

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^32 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 Wikipedia.

III-D-3 - Classe Random

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

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

Au moment de créer l'objet Random, c'est-à-dire quand la classe est instanciée, on met à jour ses attributs.

Cette méthode permet ainsi d'initialiser les attributs m, a, et c de la classe, puis le générateur de nombres aléatoires en mettant à jour x à l'aide de l'argument s.

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
)

Les valeurs par défaut de m, a et c sont également celles utilisées par le générateur d'UNIX.

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

A l'instar de la fonction random() du module random, elle permet de générer un nombre aléatoire compris entre 0 et 1.

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

Les attributs m, a, c et x sont ainsi accessibles depuis toutes les méthodes de la classe et x peut également être mis à jour par ces fonctions.

La classe Random contient également une méthode init_rand() pour initialiser le générateur de nombres aléatoires une fois l'objet créé, et une méthode rand_int(lower,upper) permettant de générer un entier au hasard. Libre à chacun d'ajouter d'autres méthodes à cette classe.

III-D-4 - Test du générateur à l'aide d'une simulation de Monte-Carlo

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 Wikipedia.

III-D-4-a - Fonction Test

Soit un point M de coordonnées (x,y), où et . On tire aléatoirement les valeurs de x et y avec notre méthode rand(). Le point M appartient au disque de centre (0,0) de rayon 1 si et seulement si . La probabilité que le point M appartienne au disque est donc Pi/4.

Donnons maintenant le déroulé de la fonction test :

- 1. Création de l'objet Random et initialisation du générateur de nombres aléatoires ;

- 2. Pour i=0 jusqu'à n-1 faire :

- --- 2.1. Tirage aléatoire des coordonnées (x,y) du point M et ajout du point à la liste ;

- --- 2.2. Incrémentation du compteur si le point est dans le cercle ;

- 3. Renvoi de la valeur de (cpt / n) * 4 comme approximation de Pi et de la liste des points.

On obtient ainsi le code :

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)

On crée l'objet Random au début ce qui initialise aussi le générateur de nombres aléatoires.

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é :

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 :

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

Puis :

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

Et enfin :

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.

IV - Conclusion

Nous voici donc arrivé au terme de notre article. Nous avons donc d'abord appris à utiliser les principales fonctions du module random. Dans un second temps et à l'aide de ces fonctions, nous avons également su réaliser différentes simulations (jeu de pile ou face, échantillonnage aléatoire, etc.). Enfin, nous avons proposé au lecteur de créer et de tester son propre générateur de nombres pseudo-aléatoires à l'aide d'une méthode de Monte-Carlo.

Ces algorithmes restent pour le moment plus simples à mettre en œuvre et plus efficaces que les générateurs de « vrais » nombres aléatoires. Cependant, il faut rappeler que dans certaines circonstances, l'utilisation de nombres pseudo-aléatoires à la place de « vrais » nombres aléatoires peut totalement compromettre l'objectif recherché. C'est le cas par exemple en cryptologie où les clefs de chiffrement doivent être parfaitement aléatoires pour garantir une sécurité maximale. Il existe dans ce cas des librairies et des fonctions plus adaptées, comme par exemple le module secrets, ou encore la fonction urandom() du module os.

VI - Téléchargement

Le fichier compressé contenant le module Python pour effectuer les différents tests.

VII - Remerciements

Je tiens à remercier f-leb et Laurent Ott pour leurs différentes sugggestions, ainsi que ... pour sa relecture.

Sources

- Page Wikipédia sur les générateurs de nombres pseudo-aléatoires

- Page sur les graines aléatoires

- Documentation sur le module random

- Page sur les échantillons en statistiques

- Générateur congruentiel linéaire

- Congruence sur les entiers

- Méthode de Monte-Carlo