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 kitxmlcodeinlinelatexdvp\overline{X}_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvpXfinkitxmlcodeinlinelatexdvp est une variable aléatoire dont l'univers image est kitxmlcodeinlinelatexdvpX(Ω) = \left\{ x_{1},x_{2},\cdots ,x_{n} \right\}finkitxmlcodeinlinelatexdvp et que l'on connait kitxmlcodeinlinelatexdvpp_{1} = P(X=x_{1}),p_{2} = P(X=x_{2}),\cdots , p_{i} = P(X=x_{i}),\cdots , p_{n} = P(X=x_{n})finkitxmlcodeinlinelatexdvp, avec kitxmlcodeinlinelatexdvpp_{1} + p_{2} + \cdots + p_{n} = 1finkitxmlcodeinlinelatexdvp, alors l'ensemble des couples kitxmlcodeinlinelatexdvp\left( x_{i} , p_{i} \right)finkitxmlcodeinlinelatexdvp est la distribution de probabilité de kitxmlcodeinlinelatexdvpXfinkitxmlcodeinlinelatexdvp, aussi appelée loi de probabilité de kitxmlcodeinlinelatexdvpXfinkitxmlcodeinlinelatexdvp.
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 :
kitxmlcodelatexdvpf(x) = \frac{1}{\sigma\sqrt{2\pi}} \ e^{\textstyle -\frac{1}{2} {\left(\frac{x - \mu}{\sigma}\right)}^2}finkitxmlcodelatexdvpParmi 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 kitxmlcodeinlinelatexdvpX_{1},X_{2},⋯,X_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp{S}_{n}=X_{1}+X_{2}+⋯+X_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}=S_{n}/nfinkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}finkitxmlcodeinlinelatexdvp appartiennent à l'intervalle kitxmlcodeinlinelatexdvp\left[ \mu - 2\sigma/\sqrt{n}, \mu + 2\sigma/\sqrt{n} \right]finkitxmlcodeinlinelatexdvp.
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 kitxmlcodeinlinelatexdvpX_{1},X_{2},⋯,X_{n}finkitxmlcodeinlinelatexdvp 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 :
Par conséquent, d'après le théorème central limite, la variable kitxmlcodeinlinelatexdvp\overline{X}_{n}=(X_{1}+X_{2}+⋯+X_{n})/nfinkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}finkitxmlcodeinlinelatexdvp, 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 kitxmlcodeinlinelatexdvp{S}_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvpx_{i} = x_{0} + ih , i \in \left\{ 1 , … , n \right\}, \forall h \neq 0finkitxmlcodeinlinelatexdvp, on peut alors estimer à l'aide de la fonction précédente la proportion de valeurs comprises dans l'intervalle kitxmlcodeinlinelatexdvp\left[x_{i}-h/2, x_{i}+h/2\right]finkitxmlcodeinlinelatexdvp et donc égale à kitxmlcodeinlinelatexdvpx_{i}finkitxmlcodeinlinelatexdvp :
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 :
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 kitxmlcodeinlinelatexdvpX_{1},X_{2},⋯,X_{n}finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}=S_{n}/nfinkitxmlcodeinlinelatexdvp 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 :
- Initialisation des paramètres k et n, ainsi que des valeurs mini et maxi en abscisse ;
- Création de la liste des valeurs en abscisse ;
- Génération des k séquences d'entiers aléatoires et des k sommes d'entiers ;
- 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 :
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 kitxmlcodeinlinelatexdvpX_{1},X_{2},⋯,X_{n}finkitxmlcodeinlinelatexdvp distribuée chacune suivant la même loi uniforme de support kitxmlcodeinlinelatexdvp\left\{ 1, 2, ..., 100 \right\}finkitxmlcodeinlinelatexdvp, d'espérance kitxmlcodeinlinelatexdvp\mu = \left( 100 + 1 \right)/2finkitxmlcodeinlinelatexdvp et de variance kitxmlcodeinlinelatexdvpv=\left( 100^{2} - 1 \right)/12finkitxmlcodeinlinelatexdvp.
Par conséquent, d'après le théorème central limite, on s'attend à ce que la variable kitxmlcodeinlinelatexdvp\overline{X}_{n}=(X_{1}+X_{2}+⋯+X_{n})/nfinkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\left[ \mu - 3\sigma/\sqrt{n}, \mu + 3\sigma/\sqrt{n} \right]finkitxmlcodeinlinelatexdvp.
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 kitxmlcodeinlinelatexdvp(\overline{X}_{n})finkitxmlcodeinlinelatexdvp 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 kitxmlcodeinlinelatexdvp\overline{X}_{n}=(X_{1}+X_{2}+⋯+X_{n})/nfinkitxmlcodeinlinelatexdvp 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éfichier de test 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