Table des matières PDFPython

Principe des méthodes de Monte-Carlo

1. Introduction

D'une manière générale, les méthodes de Monte-Carlo utilisent des nombres aléatoires pour simuler des phénomènes comportant une ou plusieurs variables aléatoires. Le nom provient du célèbre casino de Monte-Carlo.

On considère une simulation de Monte-Carlo élémentaire, visant à évaluer l'espérance (et la variance) d'une variable aléatoire en générant un grand nombre d'échantillons qui suivent la même loi de probabilité que la variable aléatoire.

2. Variable aléatoire discrète

2.a. Espérance et variance

On considère une variable aléatoire X pouvant prendre les M valeurs xk avec k=0,..M-1. La probabilité de l'évènement xk est notée pk. La donnée de ces probabilités constitue la loi de la variable aléatoire, appelée aussi distribution des probabilités. La somme des probabilités doit être égale à 1 :

L'espérance de la variable aléatoire est :

Remarque : en physique statistique, l'espérance d'une grandeur physique aléatoire est appelée moyenne ou moyenne statistique de cette grandeur.

La variance est l'espérance du carré de l'écart entre la variable et son espérance :

L'écart type est la racine carrée de la variance :

2.b. Simulation de Monte-Carlo

Principe

On utilise un générateur de nombres pseudo-aléatoires pour générer des échantillons xi de la variable aléatoire X. Les valeurs de ces échantillons sont dans l'ensemble des valeurs xk.

Soit N le nombre d'échantillons générés. L'espérance est évaluée en calculant la moyenne empirique définie par :

Cette moyenne est généralement calculée pour des valeurs de N croissantes, ce qui nécessite le stockage de la somme des xi.

La variance empirique est définie par :

L'estimation de la variance nécessite donc de stocker la somme des carrés des échantillons.

La moyenne yM est elle-même une variable aléatoire : pour N fixé, sa valeur varie aléatoirement d'un tirage à l'autre. La variance de yN est d'autant plus faible que N est grand. Pour calculer cette variance, on utilise les deux propriétés suivantes valables pour deux variables aléatoires indépendantes :

On obtient ainsi la variance de la moyenne :

Il s'en suit que l'écart-type sur la moyenne varie comme l'inverse de la racine carrée :

Cet écart-type peut être évalué en utilisant l'évaluation vN de la variance de X.

Le théorème de la limite centrale ([1]) établit que la moyenne empirique suit (dans la limite N infini) une distribution continue gausienne (voir plus loin). Un intervalle de confiance à 95 pour cent est donné par :

La probabilité pour que l'espérance cherchée se trouve dans cet intervalle est de 0,95.

Exemple

On considère un tirage de dé. Les valeurs possibles sont 1,2,3,4,5,6 avec chacun la probabilité 1/6. La fonction random.randint(1,6) permet d'obtenir les échantillons. Dans ce cas élémentaire, on connait la valeur exacte de l'espérance :

et de la variance :

La fonction suivante effectue N tirages, calcule la moyenne et la variance empiriques, et le demi-intervalle de confiance à 95 pour cent :

import random
import math

def tirages(N):
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        x = random.randint(1,6)
        somme += x
        somme2 += x*x
    moyenne = somme*1.0/N
    variance = somme2*1.0/N-moyenne*moyenne
    ecart = 1.96*math.sqrt(variance*1.0/N)
    return (moyenne,variance,ecart)
                       

Voici un exemple avec 1000 tirages :

t = tirages(1000)
                       
print(t)
--> (3.431, 2.889239, 0.10535321799736352)

Voici une deuxième série de tirages :

t = tirages(1000)
                       
print(t)
--> (3.462, 2.9685559999999995, 0.10678953473819426)

Pour réduire l'écart-type de la moyenne, on doit augmenter le nombre de tirages :

t = tirages(10000)
                       
print(t)
--> (3.4791, 2.8991631900000012, 0.033372781290602685)

Cette dernière simulation donne donc, avec un intervalle de confiance à 95 pour cent :

3. Variable aléatoire continue

3.a. Densité de probabilité

Soit X une variable aléatoire continue, qui peut prendre toute valeur réelle dans l'intervalle [a,b]. On définit une fonction p(x) appelée densité de probabilité de la variable aléatoire, telle que la probabilité d'obtenir une valeur dans l'intervalle [x1,x2] soit :

La densité de probabilité doit vérifier (si possible) la condition suivante, appelée condition de normalisation :

On introduit aussi la fonction de répartition F(x) qui donne la probabilité d'obtenir une valeur inférieure ou égale à x :

En dérivant la fonction de répartition, on obtient :

La probabilité d'obtenir une valeur x donnée dans l'intervalle [a,b] est nulle. En pratique, on s'intéresse à la probabilité d'obtenir une valeur entre x et x+δx, égale approximativement à :

densite.svgFigure pleine page

L'espérance d'une variable aléatoire f(X) est définie par :

La variance se définit comme pour une variable discrète.

Une densité de probabilité uniforme est constante sur l'intervalle [a,b]. Si les bornes de cet intervalle ne sont pas à l'infini, la densité de probabilité est obtenue avec la condition de normalisation :

L'espérance de x est alors :

La distribution gaussienne (ou normale) est définie sur l'intervalle par :

Elle vérifie la condition de normalisation et on a :

3.b. Simulation de Monte-Carlo

Principe

Considérons la simulation d'une variable aléatoire continue définie sur l'intervalle [a,b]. Les nombres réels sont représentés par des nombres à virgule flottante, qui ont bien sûr une précision limitée. Si l'on note δx cette précision, l'intervalle [a,b] est divisé en M sous-intervalles égaux de largeur δx. Le tirage d'un flottant dans cet intervalle revient donc à tirer un entier compris entre 0 et M. La probabilité d'obtenir un nombre flottant x est égale à p(x)δx.

La fonction random.uniform(a,b) délivre un flottant avec une densité de probabilité uniforme sur l'intervalle [a,b]. On peut aussi utiliser la fonction random.random() qui fait la même chose sur l'intervalle [0,1[ (la valeur 1 est exclue).

La fonction random.gauss(mu,sigma) délivre un flottant avec la distribution de Gauss.

L'estimation de l'espérance et de la variance se fait comme déjà expliqué pour une variable aléatoire discrète.

Exemple

On considère comme exemple le tirage de nombres aléatoires avec la loi de Gauss. La fonction suivante fait N tirages. Elle calcule la moyenne empirique (évalutation de l'espérance), la variance empirique, et l'intervalle de confiance à 95 pour cent pour l'espérance.

def tirages_gauss(N,mu,sigma):
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        x = random.gauss(mu,sigma)
        somme += x
        somme2 += x*x
    moyenne = somme*1.0/N
    variance = somme2*1.0/N-moyenne*moyenne
    ecart = 1.96*math.sqrt(variance*1.0/N)
    return (moyenne,variance,ecart)
                     

Voici un exemple :

t = tirages_gauss(10000,1.0,0.1)
                     
print(t)
--> (0.9991884164174464, 0.00999573333003767, 0.0019595818217332164)

On obtient ainsi l'estimation de l'espérance suivante :

3.c. Échantillonnage d'une densité non uniforme

Une distribution continue de densité non uniforme peut être échantillonnée en inversant la fonction de répartition, définie par :

On suppose que la densité de probabilité est strictement positive. La fonction de répartition est alors continue et strictement croissante. De plus F(a)=0 et F(b)=1. Il s'en suit qu'elle admet une fonction inverse F-1 définie sur l'intervalle [0,1] et à valeurs dans [a,b]. Considérons alors une variable aléatoire U de densité uniforme sur l'intervalle [0,1] et soit la variable X définie par :

La probabilité d'obtenir x compris entre x1 et x2 est :

ce qui montre que X a la densité p. Si la fonction de répartition est inversible analytiquement, on peut obtenir l'échantillonnage à partir d'un échantillonnage uniforme sur l'intervalle [0,1].

Considérons par exemple une densité de probabilité proportionnelle à x sur l'intervalle [0,1] :

En écrivant la condition de normalisation on obtient :

La fonction de répartition est :

L'inversion s'écrit :

u est tiré aléatoirement avec une densité uniforme sur l'intervalle [0,1].

Pour tester l'échantillonnage, on fait un histogramme avec plusieurs milliers de tirages :

import numpy.random
from matplotlib.pyplot import *

N = 100000
x = numpy.sqrt(numpy.random.random_sample(N))
hist(x,100)

                 
figA.svgFigure pleine page

Si l'inverse de la fonction de répartition n'est pas disponible, on peut la calculer numériquement et la stocker dans une table (sous forme discrète).

Pour la génération de nombres suivant une loi discrète, voir Échantillonnage des distributions de probabilité discrètes.

Références
[1]  B. Jourdain,  Probabilités et statistique,  (Ellipses, 2009)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.