Principe des méthodes de Monte-Carlo

Principe

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

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

$$y_N=moy(X,N)=\frac{1}{N}\sum _{i=0}^{N-1}{x}_i\text{(8)}$$

L'espérance de la moyenne empirique est égale à l'espérance de X.

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

La variance empirique est définie par :

$$v_N=var(X,N)=\frac{1}{N-1}\sum _{i=0}^{N-1}{(x_i-y_N)}^2\text{(9)}$$

L'espérance de la variance empirique est égale à la variance de X.

En pratique, on préfère utiliser la forme développée suivante :

$$v_N=\frac{1}{N-1}\sum _{i=0}^{N-1}{x}_i^2-\frac{N}{N-1}{y}_N^2\text{(10)}$$

Dans une simulation de Monte-Carlo, le nombre N de tirage est très grand (au moins 1000). On peut donc confondre N-1 et N, ce qui donne l'expression approchée suivante de la variance empirique :

$$v_N\approx \frac{1}{N}\sum _{i=0}^{N-1}{x}_i^2-y_N^2\text{(11)}$$

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

La moyenne y_M est elle-même une variable aléatoire : pour N fixé, la répétition de N tirages de la variable X donne des valeurs aléatoires de la moyenne. La variance de y_N 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 :

$$\begin{array}{ccc}\hfill & var(x_1+x_2)=var\left(x_1\right)+var\left(x_2\right)\hfill & \hfill \text{(12)}\\ \hfill & var\left(ax\right)=a^{2}var\left(x\right)\hfill & \hfill \text{(13)}\end{array}$$

On obtient ainsi la variance de la moyenne :

$$var\left(y_N\right)=\frac{var\left(X\right)}{N}\text{(14)}$$

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

$$\Delta \left(y_N\right)=\frac{\sqrt{var\left(X\right)}}{\sqrt{N}}\text{(15)}$$

Cet écart-type est évalué avec la variance empirique v_N.

Le théorème de la limite centrale établit que la moyenne empirique suit (si N est assez grand) une distribution continue gaussienne (voir plus loin). Un intervalle de confiance à 95 pour cent est donné par :

$$[y_N-\frac{1,96\sqrt{v_N}}{\sqrt{N}},y_N+\frac{1,96\sqrt{v_N}}{\sqrt{N}}]\text{(16)}$$

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

Exemple

On considère des tirages d'un 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 connaît la valeur exacte de l'espérance :

$$E\left(X\right)=\frac{1}{6}(1+2+3+4+5+6)=\frac{7}{2}=3,5\text{(17)}$$

et de la variance :

$$var\left(X\right)=\frac{1}{6}({(1-7/2)}^2+{(2-7/2)}^2+\cdots )\simeq 2,9167\text{(18)}$$

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.526, 2.979324000000002, 0.10698304107848126)

Voici une deuxième série de tirages :

t = tirages(1000)
                       

print(t)
--> (3.418, 2.945276, 0.10636997829086928)

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

t = tirages(10000)
                       

print(t)
--> (3.5017, 2.8911971100000002, 0.03332690027256661)

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

$$E\left(X\right)=3\mathrm{.}50\pm 0\mathrm{.}03\text{(19)}$$

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 nombre à virgule flottante dans cet intervalle revient donc à tirer un entier compris entre 0 et M. La probabilité d'obtenir un nombre 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, la variance empirique, et l'intervalle de confiance à 95 pour cent pour l'espérance.

La fonction génère aussi un histogramme des échantillons. Pour générer un histogramme, il faut choisir un intervalle [a,b] et des classes de valeur dans cet intervalle. Supposons que ces classes soient les n_b sous intervalles de largeur h=(b-a)/n_b. Notons H_j le nombre de tirages donnant une valeur dans l'intervalle [a+jh,a+(j+1)h] avec j variant de 0 à n_b-1. Lorsque les valeurs de la variable aléatoire ne sont pas dans un intervalle borné, il faut choisir des valeurs de a et b. Une première approche consiste à générer tous les échantillons en les stockant dans un tableau puis à choisir le minimum de le maximum des valeurs des échantillons pour a et b. Il est préférable de fixer la largeur des classes (h) et donc de calculer n_b en conséquence. Cette approche a deux inconvénients : elle nécessite le stockage des échantillons; elle ne permet pas de générer un histogramme qui se met à jour au fûr et à mesure de la génération des échantillons. Nous préférons calculer l'histogramme en fixant a priori les valeurs de a et b. L'inconvénient de cette approche est le risque d'avoir des valeurs d'échantillons en dehors de cet intervalle. L'intervalle doit donc être choisi assez grand pour que ce risque soit négligeable.

Le nombre de tirages H_j dans une classe est finalement divisé par le nombre de tirages total (comptés dans l'histogramme) afin d'obtenir une évaluation de la probabilité d'obtenir un tirage dans cette classe.

import numpy
                     
def tirages_gauss(N,mu,sigma,a,b,nb):
    somme = 0.0
    somme2 = 0.0
    h = (b-a)/nb
    H = numpy.zeros(nb)
    n = 0
    for i in range(N):
        x = random.gauss(mu,sigma)
        somme += x
        somme2 += x*x
        j = int((x-a)/h)
        if j>=0 and j<nb:
            H[j] += 1
            n += 1
    moyenne = somme*1.0/N
    variance = somme2*1.0/N-moyenne*moyenne
    ecart = 1.96*math.sqrt(variance*1.0/N)
    H=H/n
    x=numpy.linspace(a,b,nb)
    return (moyenne,variance,ecart,x,H)
                     

Voici un exemple :

mu=1.0
sigma=0.1
a=0.5
b=1.5
nb=50
(m,v,e,x,H) = tirages_gauss(10000,mu,sigma,a,b,nb)
                     

print((m,v,e))
--> (1.0006659135758849, 0.010047662063773277, 0.001964665329876603)

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

$$E\left(X\right)=1\mathrm{.}001\pm 0\mathrm{.}002\text{(32)}$$

Voici une représentation graphique de l'histogramme :

from matplotlib.pyplot import *
figure()
plot(x,H,"o")
grid()
xlabel("x")
ylabel("P")

                     

Figure pleine page

Un histogramme d'une variable aléatoire continue permet d'évaluer la densité. Pour cela, il faut diviser chaque probabilité par la largeur du sous-intervalle correspondant. Dans le cas présent, on peut faire une comparaison avec la densité de probabilité :

h=(b-a)/nb
H=H/h
figure()
plot(x,H,"o")
x=numpy.linspace(a,b,1000)
p=1/(numpy.sqrt(2*numpy.pi)*sigma)*numpy.exp(-(x-mu)**2/(2*sigma**2))
plot(x,p,"k-")
grid()
xlabel("x")
ylabel("p")
                     

Figure pleine page