Calcul d'intégrales : méthode de Monte-Carlo

3.a. Principe

On reprend l'exemple de l'intégrale d'une fonction à une variable sur l'intervalle [a,b].

Soit X une variable aléatoire réelle ayant une densité de probabilité uniforme sur l'intervalle [a,b] :

$$p\left(x\right)=\frac{1}{b-a}\text{(6)}$$

Soient N échantillons x_i de cette variables aléatoire. L'évaluation de l'intégrale par la méthode de Monte-Carlo consiste (dans sa forme élémentaire) à calculer la somme suivante :

$$S_N=(b-a)\frac{1}{N}\sum _{i=0}^{N-1}f\left(x_i\right)\text{(7)}$$

Il s'agit de la même formule que celle de la méthode des rectangles, avec des points répartis aléatoirement sur l'intervalle.

Cette somme est une évaluation de l'espérance de f(X), multipliée par (b-a).

Pour obtenir une expression plus générale, considérons une variable aléatoire réelle X sur l'intervalle [a,b] avec une densité de probabilité p(x).

Pour une fonction g, l'espérance de g(X) est :

$$E\left(g\left(X\right)\right)={\int }_a^{b}p\left(x\right)g\left(x\right)dx\text{(8)}$$

En posant f(x)=p(x)g(x), on peut écrire :

$$E(f/p)={\int }_a^{b}f\left(x\right)dx\text{(9)}$$

Une approximation de l'intégrale est donc obtenue en évaluant l'espérance de f/p avec N échantillons x_i tirés dans l'intervalle [a,b] avec une densité de probabilité p, soit :

$$S_N=moy(f/p,N)=\frac{1}{N}\sum _{i=0}^{N-1}\frac{f\left(x_i\right)}{p\left(x_i\right)}\text{(10)}$$

Dans le cas d'une densité uniforme, on retrouve l'expression $\mathrm{(7)}$.

Notons σ² la variance de f(X)/p(X), définie par :

$${\sigma }^2=E\left({(f/p)}^2\right)-E{(f/p)}^2\text{(11)}$$

La variance de la moyenne empirique est :

$$var\left(S_N\right)=\frac{{\sigma }^2}{N}\text{(12)}$$

Il s'en suit que l'écart-type se réduit comme l'inverse de la racine carré de N lorsqu'on augmente le nombre d'échantillons. Ce résultat reste valable pour une intégrale double, triple, etc. L'évolution de la variance de la somme avec N est indépendante de la dimension d de l'espace. Pour une intégrale simple, cette convergence en $1/\sqrt{N}$ est plus lente que la convergence en 1/N de la méthode des rectangles, qui est la méthode de quadrature la plus lente. La méthode de Monte-Carlo peut cependant devenir avantageuse pour les intégrales triples ou en dimension d supérieure à 3, lorsque la précision souhaitée est faible.

Pour évaluer la précision du résultat, il faut évaluer la variance σ². Cela peut se faire avec le calcul de la moyenne, en calculant la variance empirique :

$$v_N=var(f/p,N)=\frac{1}{N}\sum _{i=0}^{N-1}(\frac{f\left(x_i\right)}{p\left(x_i\right)}{)}^2-S_N^2\text{(13)}$$

En divisant par N, on obtient la variance de la somme. L'intervalle de confiance à 95 pour cent est :

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

L'intégrale considérée a une probabilité 0,95 de se trouver dans cet intervalle.

La dépendance de la variance de la somme en fonction de N est toujours $1/\sqrt{N}$, quelle que soit la dimension de l'espace et la fonction intégrée. Pour réduire la variance, il faut réduire la variance σ² de f/p. L'amélioration de la convergence consiste donc à choisir une densité de probabilité p(x) qui minimise la variance et qui ne soit pas trop difficile à échantillonner. Supposons que $f\left(x\right)\ge 0$ et considérons la densité de probabilité particulière suivante :

$$p_o\left(x\right)=\frac{f\left(x\right)}{{\int }_a^{b}f\left(x\right)dx}\text{(15)}$$

La variance est dans ce cas :

$${\sigma }^2=var\left({\int }_a^{b}f\left(x\right)dx\right)=0\text{(16)}$$

La variance est nulle, mais cette densité de probabilité p_o(x) n'est pas réalisable car il faudrait pour cela connaître l'intégrale que l'on cherche. Pour réduire la variance, on cherche à définir une densité de probabilité qui s'approche au mieux de la forme de la fonction f, tout en étant simple et peu coûteuse à échantillonner. C'est ce qu'on appelle l'échantillonnage préférentiel (importance sampling). La densité de probabilité doit être grande là où la fonction f a une valeur absolue grande.

3.b. Exemple à une dimension

import numpy.random
import random
import math
import numpy
                

La fonction suivante effectue l'intégration de Monte-Carlo d'une fonction sur un intervalle [a,b], avec une densité de probabilité p(x) uniforme. Elle renvoie l'estimation de l'intégrale et le demi-intervalle de confiance à 95 pour cent.

def integration1d(fonction,a,b,N):
    x = a+(b-a)*numpy.random.random_sample(N)
    p = 1.0/(b-a)
    f = fonction(x)
    moyenne = f.sum()/(N*p)
    g = f*f
    variance = g.sum()*1.0/(N*p*p)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                 

Voici un exemple, avec l'intégrale déjà considérée plus haut :

def f(x):
    return numpy.sin(x)
(integrale,intervalle) = integration1d(f,0,math.pi,1000)
                 

print((integrale,intervalle))
--> (1.9776925716412701, 0.06172322830900628)

Pour le même nombre de points calculés, la méthode des rectangles est beaucoup plus précise. Pour réduire la variance, on peut chercher à augmenter la densité de probabilité autour de π/2, où la fonction est maximale. Bien sûr, ce type d'amélioration suppose que l'on ait des informations sur la fonction à intégrer. Une densité constante par morceaux est un bon choix, car elle peut être échantillonnée simplement. Pour obtenir une distribution non uniforme simple, on peut diviser l'intervalle en trois parts égales, et attribuer une probabilité double à l'intervalle central :

Figure pleine page

Pour échantillonner cette distribution, il suffit de tirer un nombre aléatoire p parmi 1,2,3,4 avec une égale probabilité. Si p=1, on tire x dans le premier intervalle avec une densité uniforme. Si p=2 ou 3, on tire x dans le second intervalle. Si p=4, on tire x dans le troisième intervalle. Pour chaque tirage, on doit diviser la valeur de f(x) par la densité de probabilité. Les densités de probabilité des trois intervalles sont respectivement égales à 1/4, 1/2 et 1/4 divisés par la largeur de chaque intervalle.

def integration1d_importance(fonction,a,b,N):
    delta = (b-a)*1.0/3
    c1 = a+delta
    c2 = a+2*delta
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        p = random.randint(1,4)
        if p==1:
            x = random.uniform(a,c1)
            f = fonction(x)*4
        elif p==2 or p==3:
            x = random.uniform(c1,c2)
            f = fonction(x)*2
        else:
            x = random.uniform(c2,b)
            f = fonction(x)*4
        somme += f
        somme2 += f*f
    moyenne = somme*delta/N
    variance = somme2*delta*delta/N-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96) 
                   

(integrale,intervalle) = integration1d_importance(f,0,math.pi,1000)
                   

print((integrale,intervalle))
--> (1.9782931763755307, 0.04799482505984376)

On a bien une réduction de la variance pour le même nombre d'échantillons. Pour évaluer l'efficacité de la nouvelle méthode, il faudrait aussi tenir compte du surcoût dû au tirage d'un nombre aléatoire supplémentaire. Lorsque l'évaluation de f est très coûteuse (ce n'est pas le cas ici) et domine le temps de calcul, l'échantillonnage préférentiel est certainement un bon choix.

3.c. Exemple à deux dimensions

Considérons le domaine d'intégration Ω = [-1,1]x[-1,1] et le disque D de rayon 1. On cherche à calculer l'intégrale double suivante :

$${\int }_{D}f(x,y)dxdy\text{(17)}$$ Figure pleine page

Il est possible d'échantillonner avec une densité uniforme dans un disque (voir plus loin), mais une méthode plus générale consiste à échantillonner uniformément dans le domaine carré Ω. D'une manière plus générale, le domaine d'intégration D peut avoir une forme complexe, et on choisit le domaine rectangulaire Ω le plus petit (parallèle aux axes) qui contienne D. Il suffit alors de poser f(x)=0 en dehors de D.

Voici une fonction qui fait l'intégration dans un domaine rectangulaire. La densité de probabilité est l'inverse de l'aire du rectangle.

def integration2d(fonction,a,b,c,d,N):
    x = a+(b-a)*numpy.random.random_sample(N)
    y = c+(d-c)*numpy.random.random_sample(N)
    p = 1.0/((b-a)*(d-c))
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        f = fonction(x[i],y[i])
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                  

Voici par exemple le calcul de l'intégrale de la fonction qui vaut 1 sur le disque, qui est égale à π :

def f(x,y):
    if x*x+y*y <= 1:
        return 1
    else:
        return 0

(integrale,intervalle) = integration2d(f,-1,1,-1,1,10000)
                   

print((integrale,intervalle))
--> (3.156, 0.0319886411440061)

On remarque que les points tirés en dehors du cercle ne sont pas pris en compte dans la somme (puisque f prend une valeur nulle). Cela revient à échantillonner des points dans le disque avec une méthode de rejet.

Pour revenir au cas plus général d'une fonction f à intégrer sur le disque, on remarque qu'il y a un échantillonnage préférentiel évident consistant à échantillonner seulement sur le disque. Pour ce faire, il ne faut pas utiliser une méthode de rejet, qui nous ramènerait au cas précédent. Dans le cas du disque, l'inversion de la fonction de répartition a une solution analytique. Pour deux variables aléatoires (u₁,u₂) de densité uniforme sur l'intervalle [0,1], les variables polaires suivantes :

$$\begin{array}{ccc}\hfill & \theta =2\pi u_1\hfill & \hfill \text{(18)}\\ \hfill & r=R\sqrt{u_2}\hfill & \hfill \text{(19)}\end{array}$$

donnent une distribution uniforme sur le disque de rayon R. La fonction suivante calcule une intégrale sur le disque de rayon R. Les variables de la fonction sont en coordonnées polaires.

def integrationDisque(fonction,R,N):
    theta = 2*numpy.pi*numpy.random.random_sample(N)
    r = R*numpy.sqrt(numpy.random.random_sample(N))
    p = 1.0/(numpy.pi*R*R)
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        f = fonction(r[i],theta[i])
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96)
                       

Dans ce cas, le calcul de π avec f=1 n'a plus de sens puisque la densité est l'inverse de π. Considérons plutôt l'intégration d'une fonction :

def f(r,theta):
    return 1-r*r
    
(integrale,intervalle) = integrationDisque(f,1.0,10000)
                        

print((integrale,intervalle))
--> (1.5654527243187746, 0.017746448191176947)

Lorsqu'on intègre une fonction décroissante, on a intérêt à mettre en place un échantillonnage préférentiel qui favorise la zone centrale.

3.d. Exemple en dimension d

Dans les systèmes à grand nombre de degrés de libertés, on est amené à calculer des intégrales dans un espace de dimension d qui peut être très grand. On peut citer l'espace des configurations utilisé en physique statistique pour représenter la configuration d'un système comportant un grand nombre d'atomes. Les méthodes de quadrature sont très inefficaces, voire inutilisables, lorsque la dimension est élevée. Par exemple pour la méthode des rectangles, le nombre de points nécessaires pour maintenir une précision constante évolue comme N^d. Le logarithme du nombre de points est donc proportionnel à la dimension d, une loi que l'on retrouve pour toutes les méthodes de quadrature. La méthode de Monte-Carlo ne présente pas cette dépendance par rapport à la dimension. La variance de la somme calculée est :

$$var\left(S_N\right)=\frac{{\sigma }^2}{N}\text{(20)}$$

La variance σ² de la variable f/p est susceptible d'augmenter avec la dimension, mais pas selon une loi en puissance. La méthode de Monte-Carlo a donc un avantage décisif en dimension élevée.

Comme exemple, nous allons calculer, en dimension d, le volume compris entre les hypersphères de rayons R₁ et R₂. La fonction renvoie aussi l'estimation de la variance σ².

def volume(dim,R1,R2,N):
    somme = 0.0
    somme2 = 0.0
    R12 = R1*R1
    R22 = R2*R2
    p = 1.0/math.pow(2*R2,dim)
    for i in range(N):
        r2 = 0.0
        for d in range(dim):
            x = random.uniform(-R2,R2)
            r2 += x*x
        if r2>= R12 and r2<=R22:
            f = 1.0
        else:
            f = 0.0
        somme += f
        somme2 += f*f
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,variance,math.sqrt(variance/N)*1.96)
                      

On commence par la sphère de rayon 1 en dimension 3, dont le volume est 4π/3 :

(m,v,e) = volume(3,0,1,10000)
                      

print((m,v,e))
--> (4.1168, 15.98635776, 0.07836656938441033)

Voyons la sphère en dimension 5 :

(m,v,e) = volume(5,0,1,10000)
                      

print((m,v,e))
--> (5.4112, 143.87731456, 0.2350997855408839)

Comme on le voit sur cet exemple, la variance dépend de la dimension. Augmentons le nombre d'échantillons pour obtenir un écart comparable à celui de la dimension 3. Comme l'écart évolue comme $1/\sqrt{N}$, il faut augmenter N d'un facteur 90 environ :

(m,v,e) = volume(5,0,1,10000*90)
                      

print((m,v,e))
--> (5.253582222222223, 140.51450494546174, 0.024490372761522387)

Avec la méthode des rectangles, il aurait fallu augmenter le nombre de points d'un facteur N² pour maintenir la précision constante. Même si la variance augmente avec la dimension, la méthode de Monte-Carlo est incomparablement plus efficace en dimension élevée.

Voyons le volume d'une coque d'épaisseur un dixième du rayon :

(m,v,e) = volume(5,0.9,1,10**5)
                       

print((m,v,e))
--> (2.1328, 63.70076416, 0.04946846021426743)

Voici les mêmes calculs en dimension 7 (beaucoup plus longs) :

(m,v,e) = volume(7,0,1,10**7)
                       

print((m,v,e))
--> (4.7254144, 582.5235019482726, 0.014959352543089837)

(m,v,e) = volume(7,0.9,1,10**7)
                       

print((m,v,e))
--> (2.4641664, 309.34118315311105, 0.010901215937687829)

On voit apparaître une propriété de la couche sphérique lorsque la dimension augmente : son volume se rapproche du volume de la sphère complète.