Approcher une distribution normale à l'aide d'une méthode de Monte-Carlo et de la bibliothèque Matplotlib

Denis Hulo

Objectif : approcher une distribution de Gauss (ou de Laplace-Gauss) à l'aide d'une méthode de Monte-Carlo et de la bibliothèque Matplotlib.

Niveau requis : confirmé

Commentez cet article : Commentez

- Titre : Approcher une distribution normale à l'aide d'une méthode de Monte-Carlo et de la bibliothèque Matplotlib

- Auteur : Denis Hulo

- Parution : 18 janvier 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. Loi normale

- III. Théorème central limite

- IV. Simulation d'une variable de Gauss

- V. Implémentation en Python

- V-A. Fonctions utilisées

- V-A-1. Distribution de Gauss

- V-A-2. Générateur d'entiers aléatoires

- V-A-3. Fonction de lissage de la courbe

- V-B. Approximation d'une distribution de Gauss à l'aide d'une méthode de Monte-Carlo

- V-B-1. Script Python

- V-B-2. Tests du code

- VI. Conclusion

- VII. Téléchargement

- VIII. Remerciements

I - Introduction

Après avoir présenté le théorème central limite, on va montrer comment simuler à l'aide d'une méthode de Monte-Carlo, une variable aléatoire suivant une loi normale. Puis, on va implémenter ce procédé aléatoire en Python afin de vérifier visuellement sur un graphique que la variable suit approximativement une loi de Laplace-Gauss quand n est assez grand.

Un exemple simple de variable aléatoire est donné par le résultat de plusieurs lancers d’un dé à six faces, pour lequel les valeurs possibles sont 1, 2, 3, 4, 5 ou 6. Un autre exemple est la variable aléatoire qui associe au résultat des jets de n dés la somme de leurs valeurs.

Distribution de probabilité :

Si est une variable aléatoire dont l'univers image est et que l'on connait , avec , alors l'ensemble des couples est la distribution de probabilité de , aussi appelée loi de probabilité de .

II - Loi normale

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 - Théorème central limite

Soit une suite de variables aléatoires réelles définies sur le même espace de probabilité, indépendantes et identiquement distribuées suivant la même loi D. Supposons que l'espérance μ et l'écart-type σ de D existent et soient finis avec σ ≠ 0.

Alors la variable somme a pour espérance nμ et pour écart-type σ√n

De plus, quand n est assez grand, la loi normale N(nμ, nσ^2) est une bonne approximation de la loi de Sn.

Finalement, on peut également en déduire que la variable suit approximativement une loi normale d'espérance μ et d'écart-type σ/√n quand n est suffisamment grand.

Pour être précis, dans ce cas près de 95.45 % des valeurs de appartiennent à l'intervalle .

IV - Simulation d'une variable de Gauss

Soit n tirages avec remise d'entiers consécutifs compris entre 1 et 100, chaque entier ayant la même probabilité d'être choisi. On va utiliser un générateur de nombres pseudo-aléatoires pour simuler ces tirages.

Les variables aléatoires qui représentent chaque résultat des n tirages suivent donc toutes la même loi uniforme. De plus, les tirages étant successifs, les résultats n'ont pas de lien de dépendance entre eux et ainsi les variables aléatoires sont indépendantes.

Si maintenant on réitère k fois ces n tirages d'entiers consécutifs, avec k et n suffisamment grands, on peut donc s'attendre à avoir une répartition à peu près uniforme des effectifs d'entiers pour chaque variable :

Tableau de répartition pour 100 000 séquences de 100 entiers aléatoires

Par conséquent, d'après le théorème central limite, la variable devrait suivre approximativement une loi normale si n est suffisamment grand.

Cette méthode va nous permettre également de calculer la moyenne et l'écart-type des valeurs observées de , afin ensuite de les comparer aux paramètres attendus μ et σ/√n de la loi normale.

Dans le cas de n tirages de k entiers choisis entre 0 et 1 (ex. : jeu de pile ou face), la variable somme représente le nombre de succès (ou de face) dans une suite d'épreuves de Bernoulli (0 : pile, 1 : face) et suit une loi binomiale d'espérance n/2 et d'écart-type √n/2 qui peut également être approchée par une loi normale.

V - Implémentation en Python

V-A - Fonctions utilisées

V-A-1 - Distribution de Gauss

L'objet norm du module scipy.stats représente une variable aléatoire continue suivant une loi normale qui est en fait une instance de la classe rv_continuous avec ses différentes méthodes.

La fonction de densité de la loi normale norm.pdf(x, loc=μ, scale=σ) d'espérance μ et d'écart-type σ peut être invoquée comme ceci :

x = 10; μ = 50; σ = 5

val = norm.pdf(x, loc=μ, scale=σ)

Si les valeurs de la variable aléatoire sont équidistantes, c’est-à-dire , on peut alors estimer à l'aide de la fonction précédente la proportion de valeurs comprises dans l'intervalle et donc égale à  :

xi = 10; μ = 50; σ = 5

h = 0.5
val = h*norm.pdf(xi, loc=μ, scale=σ)

ou encore en utilisant la fonction de répartition de la loi normale norm.cdf(x, loc=μ, scale=σ) :

xi = 10; μ = 50; σ = 5

h = 0.5
val = norm.cdf(xi+h/2, loc=μ, scale=σ) - norm.cdf(xi-h/2, loc=μ, scale=σ)

Cette formule est presque équivalente à la précédente pour des écarts de h entre deux valeurs consécutives.

Traçons maintenant la courbe représentative de la distribution de Gauss pour des valeurs équidistantes de h et à l'aide de la bibliothèque Matplotlib :

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

# nombre d'épreuves
n = 80

# probabilité de succès
p = 0.5

# moyenne μ et écart-type σ de la loi normale
μ = n*p
σ = sqrt(n*p*(1-p))

# é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)

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

# création de la liste des valeurs en y
# y = [(norm.cdf(xi+h/2,μ,σ)-norm.cdf(xi-h/2,μ,σ)) for xi in x]

# trace la courbe correspondant à la série de valeurs
# plt.plot(x, y, color="red")
plt.bar(x, h*norm.pdf(x,μ,σ), width = h, color="white", edgecolor = 'blue')
plt.plot(x, h*norm.pdf(x,μ,σ), color="red")

# 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("Courbe représentative de la distribution normale N({0}, {1})".format(round(μ,2),round(σ**2,2)))

# affiche le graphique
plt.show()

Le code affiche :

Courbe de Gauss

V-A-2 - Générateur d'entiers aléatoires

La fonction randint(low, high, size) du sous-module numpy.random renvoie une séquence d'entiers aléatoires compris entre low et high-1 et où size indique le nombre d'entiers générés.

Chaque entier compris entre low et high-1 a donc la même probabilité d'être choisi, et par conséquent chaque séquence d'entiers aléatoires peut être associée à une suite de variables aléatoires indépendantes distribuée chacune suivant la même loi uniforme discrète

La fonction randint peut être invoquée comme ceci :

import numpy as np

# génération de 200 entiers choisis au hasard entre 1 et 100
seq = np.random.randint(1, 101, 200)

V-A-3 - Fonction de lissage de la courbe

La fonction savgol_filter du module scipy.signal permet de lisser une courbe et ainsi d'atténuer ce qui peut-être considéré comme une perturbation ou un bruit de mesure (en traitement du signal). Elle est basée sur l'algorithme de Savitzky-Golay qui réalise une régression pour déterminer le polynôme minimisant l'erreur au sens des moindres carrés.

Dans notre cas, on souhaite lisser les points représentant les fréquences obtenues en fonction des valeurs en x.

La fonction savgol_filter peut alors être utilisée comme ceci :

from scipy.signal import savgol_filter

# liste des fréquences obtenues
f_values = [0.04, 0.05, 0.11, 0.23, 0.27, 0.11, 0.1, 0.06, 0.03]

# taille de la fenêtre
window_length = 5

# degré du polynôme
poly_order = 2

# lissage des points :
f2_values = savgol_filter(f_values, window_length, poly_order)

Si la version de la fonction savgol_filter est relativement ancienne, l'argument window_length peut devoir être un entier positif impair.

V-B - Approximation d'une distribution de Gauss à l'aide d'une méthode de Monte-Carlo

V-B-1 - Script Python

L'objectif va être maintenant de générer un graphique en Python permettant de montrer que la variable aléatoire suit approximativement une loi de Laplace-Gauss quand n est suffisamment grand. Pour cela, on va d'abord décrire les différentes étapes conduisant à ce résultat :

1. Initialisation des paramètres k et n, ainsi que des valeurs mini et maxi en abscisse ;

2. Création de la liste des valeurs en abscisse ;

3. Génération des k séquences d'entiers aléatoires et des k sommes d'entiers ;

4. Traçage du graphique représentant la courbe de Gauss et les fréquences d'apparition de chaque valeur de x.

Reprenons maintenant plus en détail ces différents points :

1. Initialisation des paramètres k et n, ainsi que des valeurs mini et maxi en abscisse :

# nombre d'entiers par séquence
k = 1000000

# nombre de séquences ou de variables aléatoires
n = 100

...

# valeurs mini et maxi en abscisse
x_min = 1
x_max = 100

# somme mini
s_min = x_min*n

2. Création de la liste des valeurs en abscisse, et initialisation du tableau des valeurs en ordonnée :

# nombre de valeurs sur l'axe des abscisses
values_nb = (x_max - x_min)*n + 1

# écart entre deux valeurs consécutives
h = 1/n

# initialisation des listes de valeurs et de fréquences
x_values = np.arange(x_min, x_max+1/n, h) # [x_min,..., x_max]
f_values = np.array([0]*values_nb) # [0, ..., 0]

3. Génération des k séquences d'entiers aléatoires et des k sommes d'entiers aléatoires :

# initialisation du générateur aléatoire
np.random.seed()

# génération des k séquences d'entiers aléatoires et des k sommes d'entiers aléatoires (Sn)
for i in range(k):

# tirage aléatoire d'une séquence de n entiers compris entre x_min et x_max, puis réalisation de la somme des entiers obtenus
s_value = sum(np.random.randint(x_min, x_max+1, size=n))

# détermination de l'indice de la valeur dans la liste x_values
index_value = (s_value - s_min) # équivaut à : index_value = (s_value/n - x_min)*n

# incrémentation du compteur de la valeur d'indice index_value dans la liste f_values
f_values[index_value] += 1

4. Traçage du graphique représentant la courbe de Gauss et les fréquences d'apparition de chaque valeur :

# calcul de la moyenne des valeurs en tenant compte de leur poids
m = np.average(x_values, weights=f_values)

# calcul de la variance des valeurs en tenant compte de leur poids
v = np.average((x_values-m)**2, weights=f_values)

# calcul de l'écart-type
sd = math.sqrt(v)

print()
print("Moyenne observée =", m)
print("Écart-type observé =", sd)

# division des effectifs par le nombre k de séquences
f_values = f_values/k

# calcul des bornes inférieure et supérieure en abscisse
low_value = round((μ - 3*σ_n)*n)
high_value = round((μ + 3*σ_n)*n)

# lissage des points :

# taille de la fenêtre
window_length = (high_value-low_value)//2

if window_length % 2: # si impair
f2_values = savgol_filter(f_values, window_length, 4)
else:
f2_values = savgol_filter(f_values, window_length + 1, 4)

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

# génération du graphique représentant les fréquences de chaque valeur de x
plt.plot(x_values, f_values, color="grey", label='fréquences observées')

# génération de la courbe représentant les fréquences observées après lissage des points
plt.plot(x_values, f2_values, color="blue", label='lissage des points')

# génération de la courbe représentative de la distribution normale N(μ, (σ^2)/n)
plt.plot(x_values, h*norm.pdf(x_values, μ, σ_n), color="red", label='N(μ,σ' + chr(0x00B2) + '/n)')

##y_values = [(norm.cdf(x_value+h/2, μ, σ_n)-norm.cdf(x_value-h/2, μ, σ_n)) for x_value in x_values]
##plt.plot(x_values, y_values, color="red", label='N(μ,σ' + chr(0x00B2) + '/n)')

# limites sur l'axe des x : [μ - 3σ/√n, μ + 3σ/√n]
plt.xlim(μ - 3*σ_n, μ + 3*σ_n)

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

plt.title("Distribution de la variable aléatoire X̄n pour n={0}".format(n))

# affichage du graphique
plt.show()

On réalise également un lissage des points à l'aide de la fonction savgol_filter du module scipy.signal.

Code complet :

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

print("Approximation d'une distribution normale à l'aide d'une méthode de Monte-Carlo :\n")

# nombre d'entiers par séquence
k = 1000000

# nombre de séquences ou de variables aléatoires
n = 10

print("Nombre d'entiers par séquence, k =", k)
print("Nombre de séquences, n =", n)

# valeurs mini et maxi en abscisse
x_min = 1
x_max = 100

# somme mini
s_min = x_min*n

# calcul de l'espérance attendue
μ = (x_min+x_max)/2

# calcul de l'écart-type attendu
σ_n = math.sqrt(((x_max-x_min+1)**2 - 1)/(n*12))

print()
print("Espérance attendue, μ =", μ)
print("Écart-type attendu, σ/√n =", σ_n)

# nombre de valeurs sur l'axe des abscisses
values_nb = (x_max - x_min)*n + 1

# écart entre deux valeurs consécutives
h = 1/n

# initialisation des listes de valeurs et de fréquences
x_values = np.arange(x_min, x_max+1/n, h) # [x_min,..., x_max]
f_values = np.array([0]*values_nb) # [0, ..., 0]

# initialisation du générateur aléatoire
np.random.seed()

# génération des k séquences d'entiers aléatoires et des k sommes d'entiers aléatoires (Sn)
for i in range(k):

# tirage aléatoire d'une séquence de n entiers compris entre x_min et x_max, puis réalisation de la somme des entiers obtenus
s_value = sum(np.random.randint(x_min, x_max+1, size=n))

# détermination de l'indice de la valeur dans la liste x_values
index_value = (s_value - s_min) # équivaut à : index_value = (s_value/n - x_min)*n

# incrémentation du compteur de la valeur d'indice index_value dans la liste f_values
f_values[index_value] += 1

# calcul de la moyenne des valeurs en tenant compte de leur poids
m = np.average(x_values, weights=f_values)

# calcul de la variance des valeurs en tenant compte de leur poids
v = np.average((x_values-m)**2, weights=f_values)

# calcul de l'écart-type
sd = math.sqrt(v)

print()
print("Moyenne observée =", m)
print("Écart-type observé =", sd)

# division des effectifs par le nombre k de séquences
f_values = f_values/k

# calcul des bornes inférieure et supérieure en abscisse
low_value = round((μ - 3*σ_n)*n)
high_value = round((μ + 3*σ_n)*n)

# lissage des points :

# taille de la fenêtre
window_length = (high_value-low_value)//2

if window_length % 2: # si impair
f2_values = savgol_filter(f_values, window_length, 4)
else:
f2_values = savgol_filter(f_values, window_length + 1, 4)

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

# génération du graphique représentant les fréquences de chaque valeur de x
plt.plot(x_values, f_values, color="grey", label='fréquences observées')

# génération de la courbe représentant les fréquences observées après lissage des points
plt.plot(x_values, f2_values, color="blue", label='lissage des points')

# génération de la courbe représentative de la distribution normale N(μ, (σ^2)/n)
plt.plot(x_values, h*norm.pdf(x_values, μ, σ_n), color="red", label='N(μ,σ' + chr(0x00B2) + '/n)')

##y_values = [(norm.cdf(x_value+h/2, μ, σ_n)-norm.cdf(x_value-h/2, μ, σ_n)) for x_value in x_values]
##plt.plot(x_values, y_values, color="red", label='N(μ,σ' + chr(0x00B2) + '/n)')

# limites sur l'axe des x : [μ - 3σ/√n, μ + 3σ/√n]
plt.xlim(μ - 3*σ_n, μ + 3*σ_n)

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

plt.title("Distribution de la variable aléatoire X̄n pour n={0}".format(n))

# affichage du graphique
plt.show()

Les tableaux d'entiers aléatoires sont donc générés avec l'instruction np.random.randint(min_value, max_value+1, size=n).

Comme précisé plus haut, chaque séquence d'entiers aléatoires peut donc être associée à une suite de variables aléatoires indépendantes distribuée chacune suivant la même loi uniforme de support , d'espérance et de variance .

Par conséquent, d'après le théorème central limite, on s'attend à ce que la variable suive approximativement une loi normale d'espérance μ et d'écart-type σ/√n.

V-B-2 - Tests du code

On va donc maintenant tester ce code en fixant k suffisamment grand et en faisant varier n.

Test n°1

L'exécution du code affiche pour n=10 :

Approximation d'une distribution normale à l'aide d'une méthode de Monte-Carlo :

Nombre d'entiers par séquence, k = 1000000
Nombre de séquences, n = 10

Espérance attendue, μ = 50.5
Écart-type attendu, σ/√n = 9.128252844876723

Moyenne observée = 50.504093
Écart-type observé = 9.12230084394014

Génération du graphique.

On a également réalisé un lissage des points représentant les fréquences observées par valeur.

Le graphique montre que la courbe des fréquences obtenues après lissage des points (en bleu) est très proche de celle représentant la distribution de la loi normale (en rouge).

De plus, près de 99.73 % des valeurs appartiennent à l'intervalle .

Test n°2

L'exécution du code affiche pour n=30 :

Approximation d'une distribution normale à l'aide d'une méthode de Monte-Carlo :

Nombre d'entiers par séquence, k = 1000000
Nombre de séquences, n = 30

Espérance attendue, μ = 50.5
Écart-type attendu, σ/√n = 5.2701992372205435

Moyenne observée = 50.49770076666682
Écart-type observé = 5.2707908020601435

Génération du graphique.

Le graphique montre cette fois que la courbe des fréquences obtenues après lissage des points (en bleu) se superpose presque parfaitement avec celle représentant la distribution de la loi normale (en rouge). La convergence en loi de la suite vers une variable de Gauss se confirme donc bien.

Libre à chacun ensuite de vérifier cette tendance en effectuant des tests pour d'autres valeurs de n.

VI - Conclusion

Après avoir présenté le théorème central limite, on a pu montrer comment simuler une variable de Gauss à l'aide d'une méthode de Monte-Carlo. Puis, on a su implémenter ce procédé aléatoire en Python afin de vérifier sur un graphique que la variable suivait approximativement une loi normale quand n est suffisamment grand.

Chacun pourra enfin vérifier que pour n séquences d'entiers aléatoires choisis entre 1 et 100, la variable somme Sn suit à peu près une loi normale d'espérance nμ et d'écart-type σ√n conformément au théorème central limite.

VII - Téléchargement

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

VIII - Remerciements

Je tiens à remercier Malick pour le suivi de cet article, ainsi que escartefigue pour sa relecture.

Sources

- Page sur les distributions de probabilité

- Page Wikipedia sur les distributions de probabilité

- Page Wikipedia sur la loi uniforme discrète

- Loi Normale

- Théorème central limite

- Méthode de Monte-Carlo

- Générateur de nombres aléatoires

- Algorithme de Savitzky-Golay

- Problème sur le tirage cette fois sans remise de séquences de nombres entiers