1. Généralités

Les méthodes de Monte-Carlo sont des algorithmes permettant la résolution de problèmes d’optimisation, d’intégration ou d’échantillonnage. Elles reposent sur le tirage de nombres pseudo-aléatoires dont l’espérance de la distribution permet de calculer la quantité recherchée.

Ainsi, pour obtenir une approximation de $\pi$ on peut tracer aléatoirement des points de coordonnées (x, y) uniformément répartis dans un carré de côté 1. Les coordonnées de chaque point sont des réels compris entre 0 et 1. Les points appartiennent au quart de disque de centre (0, 0) de rayon r = 1 inscrit dans le carré si et seulement si $x^2+y^2\leq1$ .

La probabilité qu’un point appartienne au quart de disque est $\dfrac{\pi}{4}$.

En effet, la surface du quart de disque est $A_d=\dfrac{\pi \times r^2}{4}=\dfrac{\pi}{4}$, et celle du carré qui le contient $A_c=c^2=1$. Le rapport $\dfrac{A_d}{A_c}$ est donc $\dfrac{\pi}{4}$.

Le rapport entre le nombre de points placés dans le disque et le nombre total tend vers le rapport des aires, soit $\dfrac{\pi}{4}$, quand le nombre de points tend vers l’infini. L’approximation est d’autant plus précise que nombre de tirages est grand.

Une variante consiste à inscrire un disque complet centré en (0, 0) dans un carré de côté 2. Les coordonnées des points sont alors comprises entre -1 et 1. Le résultat est le même puisque les rapports des surfaces sont identiques :

Choix d’un algorithme

L’algorithme va consister à générer des paires de nombres aléatoires x et y compris entre 0 et 1 (ou -1 et 1 selon les cas) qui formeront les coordonnées cartésiennes des points. Il suffira ensuite de compter ceux qui remplissent la condition $x^2 + y^2\leq1$, c’est-à-dire qui sont situés à l’intérieur du cercle de centre (0, 0) et de rayon 1 puis à multiplier ce chiffre par 4 et à le diviser par le nombre total de points. Le résultat sera une approximation de $\pi$.

Si nous appelons $\pi_{mc}$ l’approximation obtenue :

$$\pi_{mc} = 4 \times \frac{ \mathrm{Nombre\ de\ points\ à\ lʼintérieur\ du\ cercle}}{\mathrm{Nombre\ de\ points\ total}}$$

Que la condition testée soit $< 1$ (strictement inférieur) ou $\leq 1$ (inférieur ou égal) est ici sans importance puisque la probabilité de tirer un point de coordonnées (1, 1) est nulle. Cela ressort naturellement du fait que l’intervalle entre 0 et 1 contient une infinité de réels. De manière générale, si l’on note X la variable aléatoire qui, au tirage d’un nombre, associe ce nombre, alors pour tout réel x de [0, 1], P(X = x) = 0.

Optimisation du traitement des données

Plus le nombre de points calculé est élevé, plus il devient indispensable d’utiliser un algorithme optimisé, sans quoi le temps de traitement peut rapidement devenir rédhibitoire.

Une approche naïve consiste à générer des points un par un, par exemple au moyen d’une boucle for.

Avec cette méthode, quel sera le temps nécessaire pour traiter un million de points ?

%%timeit
import random
acc = 0
for i in range(1_000_000):
    x = random.random()
    y = random.random()
    if (x ** 2 + y ** 2) < 1.0:
        acc += 1
pi_mc = 4.0 * acc / 1_000_000

531 ms ± 24.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

NB : %%timeit et %timeit sont des « commandes magiques intégrées » du terminal interactif iPython qui peuvent être exécutées dans le terminal ou dans l’environnement Jupyter Notebook. Elles permettent de mesurer très facilement la durée d’exécution d’une instruction ou d’une série d’instructions.

Le résultat peut être nettement amélioré si l’on utilise un compilateur à la volée (Just-in-time) comme Numba, par exemple. qui peut être installé avec le gestionnaire de paquets pip) :

$ sudo pip install numba

C’est précisément le calcul de $\pi$ par la méthode de Monte-Carlo qui est citée sur la page d’accueil de Numba comme exemple de code optimisé.

L’utilisation la plus simple de Numba consiste à ajouter un décorateur à la fonction que l’on souhaite accélérer.

# source : http://numba.pydata.org/
from numba import jit
import random

@jit(nopython=True)
def monte_carlo_pi(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

Quelle est cette fois la vitesse d’exécution de 1 000 000 de boucles ?

%timeit monte_carlo_pi(1_000_000)

17.5 ms ± 124 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Comme on le voit, Numba accélère grandement l’exécution du code, ici d’un facteur 30. Si l’on choisit un algorithme faisant appel à de nombreuses itérations, on aura tout intérêt à optimiser le code grâce à des solutions comme Numba, cffi ou cython.

Cependant Numba ne supporte que certaines instructions de Python. Le PCG64, notamment, qui est la dernière version du générateur de nombres aléatoires de NumPy que nous verrons un peu plus loin, n’est pas reconnu par Numba.

Si l’on veut gagner du temps, sans passer par un compilateur comme Numba, il faut opter pour une approche un peu différente. Sur de grands volumes de chiffres, la vectorisation sera bien plus efficace que l’itération. C’est justement ce que propose NumPy.

NumPy permet d’exécuter très rapidement des programmes comportant de grandes quantités de chiffres grâce aux opérations sur les tableaux. On ne procède plus ici par itérations (on évite d’itérer sur les tableaux) mais en utilisant des fonctions vectorisées qui s’appliquent à tous les éléments d’un tableau.

La fonction de calcul de $\pi_{mc}$ que l’on a vue plus haut peut être réécrite pour NumPy de la façon suivante :

import numpy as np

# Initialisation du générateur de nombres aléatoires
rng = np.random.default_rng() 

def mc_pi_numpy(nsamples):
    tab_xy = rng.random((nsamples, 2))
    acc = np.sum(tab_xy[:, 0]**2 + tab_xy[:, 1]**2 < 1)
    return 4.0 * acc / nsamples

Bien que le temps d’exécution soit un peu plus long qu’avec Numba, il est du même ordre de grandeur :

%timeit mc_pi_numpy(1_000_000)

27.6 ms ± 503 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Avec NumPy, comme on le voit, on se contente de créer un tableau de $n\times 2$ nombres aléatoires puis on applique à ce tableau des fonctions qui sont mises en œuvre sur l’ensemble des éléments sans qu’il soit nécessaire d’itérer.

Cependant, si le nombre de points calculés devient très grand, tous les éléments du tableau étant stockés dans la mémoire, celle-ci risque d’être assez rapidement saturée.

Répartir les données

Si l’on essaie de créer un tableau d’un milliard de points, il y a des chances que ça ne passe pas.

mc_pi_numpy(1_000_000_000) # un milliard (10⁹)

Même sur les serveurs cloud de Google, avec Colab (version non professionnelle), cette instruction provoque un crash du kernel :

Votre session a planté après avoir utilisé toute la mémoire RAM disponible. Si vous êtes intéressé par un accès à des environnements d’exécution à mémoire RAM élevée, nous vous conseillons Colab Pro.

tcmalloc: large alloc 16000000000 bytes == [...]

Cette limite peut être repoussée en allouant plus de RAM aux calculs mais elle n’en est pas moins un obstacle si l’on veut obtenir des résultats sur un grand nombre de points (au-delà d’un milliard).

Pour contourner le problème, il suffira d’ajouter une boucle qui permettra de répartir les coordonnées des points sur plusieurs tableaux. Au lieu de créer un tableau contenant 10⁹ données, on pourra parcourir 10 tableaux qui en comportent chacun 10⁸, 100 qui en contiennent 10⁷, etc.

Précision des calculs

Une autre limite à laquelle on va devoir faire face, est celle de la précision des calculs. Les nombres à virgule flottante sur 64 bits sont le dtype (Data type objects) par défaut utilisé par la fonction np.random.default_rng(). Ici, la précision n’est pas en problème en soi puisqu’il s’agit de nombres aléatoires dont on attend seulement qu’ils soient uniformément répartis.

En revanche, au moment du calcul final qui donne la valeur $\pi_{mc}$ à partir d’opérations entre trois entiers, on aura besoin d’une grande précision si l’on souhaite exploiter les résultats au-delà de quelques décimales.

On utilisera pour cela le module decimal de la bibliothèque standard, qui permet de faire des calculs sur les nombres à virgule flottante avec une précision aussi grande que souhaitée.

Supposons par exemple qu’un programme crée 26 805 949 036 points aléatoires que, sur ce nombre, 21 053 343 141 appartiennent au disque.

Il faut pouvoir interpréter ce résultat et savoir combien l’opération $4\times$$ \dfrac{21 \, 053 \,343 \, 141}{26 \, 805 \, 949 \, 036}$ va donner de décimales de $\pi$ exactes.

Par défaut, la précision n’est pas suffisante pour répondre à cette question :

print(4 * 21_053_343_141 / 26_805_949_036)

3.141592653589793

On voit simplement que l’on a obtenu au moins 15 décimales exactes (en supposant que la dernière affichée ne soit pas un arrondi), mais nous n’en savons pas beaucoup plus. C’est ici que le module decimal révèle toute son utilité.

from decimal import Decimal, getcontext
getcontext().prec = 32
getcontext().rounding = 'ROUND_DOWN'

print((4 * Decimal(21_053_343_141) / 
       Decimal(26_805_949_036)).quantize(Decimal('.'+'0'*30)))

3.141592653589793238462381742774

Dans cet exemple, on a fixé la précision à plus de 30 chiffres après la virgule et l’on constate que l’opération nous donne des valeurs de $\pi$ exactes jusqu’à la 21e décimale.
Le choix de ROUND_DOWN pour arrondir les résultats permettra par exemple d’afficher 1,666…666 au lieu de 1,666…667 comme résultat de l’opération 10 / 6. Car ce qui nous intéresse ici, c’est d’obtenir le plus grand nombre de décimales exactes, pas forcément d’afficher la valeur la plus proche au sens mathématique.
La méthode quantize() permet d’afficher le résultat avec le nombre de décimales souhaité.
La précision utilisée dans les calculs doit être un peu supérieure à la précision d’affichage, car elle doit tenir compte des erreurs accumulées lors des arrondis.

La fonction de calcul des $\pi_{mc}$ peut donc être améliorée grâce à la répartition des calculs entre plusieurs tableaux et à l’utilisation du module decimal pour le calcul final :

import numpy as np
import decimal as d
c_func = d.Context(prec=32, rounding='ROUND_DOWN')

rng = np.random.default_rng() 

def mc_pi_numpy(ntab, loops):
    d.setcontext(c_func)
    acc = 0
    for _ in range(loops):
        tab_xy = rng.random((ntab, 2))
        acc += np.sum(tab_xy[:, 0]**2 + tab_xy[:, 1]**2 < 1)
    return (4 * d.Decimal(int(acc) / 
            d.Decimal(ntab * loops))).quantize(d.Decimal('.'+'0'*30))

On peut cette fois obtenir des approximations pour un nombre de points élevé, et afficher le résultat avec toute la précision souhaitée :

# 10⁶ x 10³ soit 10⁹
pi_mc1 = mc_pi_numpy(1_000_000, 1_000) 
# 12_345_678 x 89 soit 1_098_765_342
pi_mc2 = mc_pi_numpy(1_234_567, 89)   
print(f'pi_mc1 : {pi_mc1}\npi_mc2 : {pi_mc2}')

pi_mc1 : 3.141645392000000000000000000000
pi_mc2 : 3.141692247592644113416719648137

⁂

Le générateur de nombres pseudo-aléatoires (PRNG)

Puisque la génération de nombres pseudo-aléatoires est au cœur de l’algorithme, sa vitesse d’exécution sera déterminante. Examinons quelques-unes des solutions proposées par Python pour créer des séquences aléatoires et comparons leurs durées d’exécution pour des blocs de 1 000 nombres, uniformément répartis entre 0 et 1 :

import random
import numpy as np
# algorithme Mersenne Twister de la bibliothèque standard
%timeit [random.random() for _ in range(1000)]
# NumPy - Mersenne Twister MT19937 (legacy)
%timeit np.random.rand(1000) 
# NumPy - PCG64 (nouvelle version)
rng = np.random.default_rng()
%timeit rng.random(1000)

122 µs ± 6.12 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
13.2 µs ± 32.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
9.36 µs ± 422 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

La première méthode est celle fournie par la bibliothèque standard de Python. Ce n’est évidemment pas la plus efficace pour remplir des tableaux puisqu’il faut procéder soit par itération, soit par compréhension de liste, comme dans l’exemple ci-dessus.

La seconde fait appel au Mersenne Twister MT19937 de NumPy qui a été conservé comme système hérité (legacy) par souci de compatibilité ascendante, mais qui n’est plus la façon officiellement recommandée de générer des nombres pseudo-aléatoires.

La troisième méthode est fondée sur un nouvel algorithme développé en 2014, le PCG64 de la famille des Générateurs congruentiels permutés. C’est cette version qu’il est actuellement conseillé d’utiliser, et c’est aussi la plus rapide.

Nous ne nous intéresserons ici qu’à la distribution uniforme mais toutes les distributions sont implémentées. Entre l’ancien algorithme et le nouveau, les correspondances sont indiquées sur cette page.

Le constructeur recommandé par NumPy, np.random.default_rng() est un alias de Generator(PCG64()). Initialisés avec la même graine aléatoire (seed), ils donneront les mêmes résultats :

from numpy.random import Generator, PCG64
rng1 = np.random.default_rng(seed=42)
rng2 = Generator(PCG64(seed=42))
print('rng1 :', rng1.random(5))
print('rng2 :', rng2.random(5))

rng1 : [0.77395605 0.43887844 0.85859792 0.69736803 0.09417735]
rng2 : [0.77395605 0.43887844 0.85859792 0.69736803 0.09417735]

Reproductibilité

La reproductibilité d’une séquence de nombres aléatoires est assurée par le choix d’une graine aléatoire qui peut être un entier ou une séquence d’entiers (tableau NumPy, tuple, liste…) à partir de laquelle la même suite de nombres aléatoires sera générée.

Par défaut, si l’on ne fournit pas de seed, la valeur initiale dépendra de l’entropie du système d’exploitation qui n’est en principe pas reproductible.

Voici deux exemples affichant 3 séquences de 5 nombres, l’un sans seed, l’autre avec seed :

# sans seed
for _ in range(3):
    rng1 = np.random.default_rng() 
    print(rng1.random(5))  

# avec seed
print()
for _ in range(3):
    rng2 = np.random.default_rng(seed=12345)
    print(rng2.random(5))

[0.65530925 0.1851834  0.86324729 0.39255347 0.3878601 ]
[0.47366815 0.22243056 0.32236095 0.51787121 0.79535642]
[0.88876466 0.14850372 0.12622366 0.66145    0.43773946]
    
[0.22733602 0.31675834 0.79736546 0.67625467 0.39110955]
[0.22733602 0.31675834 0.79736546 0.67625467 0.39110955]
[0.22733602 0.31675834 0.79736546 0.67625467 0.39110955]

Dans le premier cas, même réinitialisé, le générateur produit des séquences différentes à chaque fois. Les séquences ne seront pas reproductibles. Dans le second cas, initialisé avec le même seed, il produit toujours la même séquence.

Boucles et appariements

Une autre condition importante pour assurer la reproductibilité des résultats a trait à l’ordre dans lequel sont générés les nombres pseudo-aléatoires.

Chaque point étant représenté par deux coordonnées du plan x et y, il faut deux nombres aléatoires pour le caractériser.

Pour construire 6 points, la disposition suivante en deux tableaux est possible, l’un pour les abscisses, l’autre pour les ordonnées :

\[ \textnormal{X} = \left[ \begin{array} {lll} r_1 & r_2 & r_3\\ r_4 & r_5 & r_6 \end{array} \right] \hspace{.5cm} \textnormal{Y} = \left[ \begin{array} {lll} r_7 & r_8 & r_9 \\ r_{10} & r_{11} & r_{12} \end{array} \right] \]

L’ordre dans lequel seront enregistrés les éléments de la suite aléatoire $(r_1, r_2, r_3,...)$ est celui qui est naturellement suivi par NumPy si on lui donne les instructions X = rng.random((2, 3)) et Y = rng.random((2, 3)). Dans ce cas, il remplira les deux tableaux, l’un après l’autre, ligne par ligne. Les nombres de la suite aléatoire seront donc rangés dans l’ordre indiqué ci-dessus.

Si l’on choisit un tableau pour les x et un autre pour les y, les points générés auront donc pour coordonnées $P_1(r_1, r_7), P_2(r_2, r_8)... $ jusqu’à $P_{6}(r_6, r_{12}$).

Cette construction ne présente pas de problème particulier et semble même assez intuitive. Elle a néanmoins l’inconvénient d’associer les paires x et y d’une façon qui dépend du nombre d’éléments présents dans les tableaux. Si le nombre d’éléments est différent les associations seront différentes :

\[ \textnormal{X} = \left[\begin{array} {lll} r_1 & r_2\\ r_3 & r_4 \end{array} \right] \hspace{.5cm} \textnormal{Y} = \left[\begin{array} {lll} r_5 & r_6\\ r_7 & r_8 \end{array} \right] \]

Cette fois, on a $P_1(r_1, r_5)$, puis $P_2(r_2, r_6), P_3(r_3, r_7)$ et $P_4(r_4, r_8)$. Autrement dit, ce ne seront pas les mêmes points que ceux du modèle précédent puisque la taille des tableaux ayant été modifiée, les associations formant les paires (x, y) ne sont pas les mêmes.

Même si l’on a pris soin de choisir un même seed, les résultats ne seront pas reproductibles pour peu que le nombre d’éléments varie.

Pour éviter cela, on aura tout intérêt — quand la reproductibilité a une importance — à associer entre eux des nombres consécutifs de la suite pseudo-aléatoire pour former les paires (x, y). De cette façon, elles resteront constantes, indépendamment de la dimension des tableaux.

tab_xy=[]
for n in (6, 4):
    rng = np.random.default_rng(42)
    tab_xy.append(rng.random((n, 2)))
tab1, tab2 = tab_xy

Ici, quelle que soit la taille des tableaux, les paires (x, y) sont constantes :

print('tableau 1 :\n', tab1)

tableau 1 :
    [[0.77395605 0.43887844]
    [0.85859792 0.69736803]
    [0.09417735 0.97562235]
    [0.7611397  0.78606431]
    [0.12811363 0.45038594]
    [0.37079802 0.92676499]]

print('tableau 2 :\n', tab2)

tableau 2 :
    [[0.77395605 0.43887844]
    [0.85859792 0.69736803]
    [0.09417735 0.97562235]
    [0.7611397  0.78606431]]

Les 4 premiers points sont communs aux tableaux 1 et 2 puisqu’ils ont les mêmes coordonnées.

Cette stabilité se retrouve aussi dans le cas de boucles parcourant successivement plusieurs tableaux. Pour un seed déterminé, le parcours de 10 tableaux contenant les coordonnées de 10 000 points donnera le même résultat que celui de 100 tableaux contenant les coordonnées de 1 000 points.

Illustration avec cette fonction qui renvoie le nombre de « succès », c’est-à-dire le nombre de points placés dans le cercle :

import numpy as np

def in_circle(ntab, loops):
    rng = np.random.default_rng(seed=42) # on fournit un seed 
    nb_cible = 0
    for _ in range(loops):
        tab_xy = rng.random((ntab, 2))
        nb_cible += np.sum(tab_xy[:, 0]**2 + tab_xy[:, 1]**2 < 1)
    return nb_cible

# 10 tableaux de 10_000 éléments
print('Pour 10 000 x 10 :', in_circle(10_000, 10))
# 100 tableaux de 1_000 éléments
print('Pour 1 000 x 100 :', in_circle(1_000, 100))

Pour 10 000 x 10 : 78752
Pour 1 000 x 100 : 78752

La méthode `advance()`

Le générateur de nombres pseudo-aléatoires (PRNG) produit donc une séquence de nombres $r_1, r_2, r_3, r_4,...$ qui sera toujours la même pour un seed déterminé.

Il est même possible de passer directement à $r_n$ sans avoir à calculer ou afficher les éléments précédents. Grâce à la méthode advance(), une fonction peu connue du PCG64, on peut en effet sauter une partie de la séquence et appeler directement la valeur de rang n du générateur.

Dans le cas de la simulation de Monte-Carlo, cela peut-être utile si l’on connaît déjà la proportion de succès (le nombre de points dans la cible) pour un nombre total n de points donnés, cette valeur étant bien sûr associée à un seed déterminé. On peut alors commencer à créer des points aléatoires directement à partir de n+1 avec le même seed sans avoir à lancer le programme à partir de 0. C’est le principal intérêt de cette méthode.

En voici une illustration. On passe à la fonction le nombre de points placés dans la cible $(nb\_cible_1)$ sur un nombre total $nb\_total_1$, et on lui indique un autre nombre total de points $nb\_total_2$ (avec $nb\_total_2 > nb\_total_1$) pour lequel on souhaite connaître $nb\_cible_2$. La fonction renvoie la valeur de $nb\_cible_2$ après avoir traité uniquement les points situés entre $nb\_total_1$ et $nb\_total_2$ :

import numpy as np
from numpy.random import Generator, PCG64

def in_circle_plus(seed, nb_total1, nb_total2, nb_cible1):
    # on avance au rang n1*2 car chaque point compte 2 coordonnées
    rng = Generator(PCG64(seed).advance(nb_total1*2))
    # on ne calcule que les points entre n1 et n2
    tab_xy = rng.random((nb_total2 - nb_total1, 2)) 
    return nb_cible1 + np.sum(tab_xy[:, 0]**2 + tab_xy[:, 1]**2 < 1)

seed = 42                  # graine aléatoire
nb_total1 = 70_000_000     # nombre total de points au départ
# nombre de points à l’intérieur du cercle pour la valeur de n1
nb_cible1 = 54_979_074
# nombre de points total pour lequel on veut connaître nb_cible2   
nb_total2 = 80_000_000     

nb_cible2 = in_circle_plus(seed, nb_total1, nb_total2, nb_cible1)
print('nb_cible2 =', nb_cible2)

nb_cible2 = 62833685

On peut s’assurer que cette valeur est bien celle que l’on aurait obtenue si l’on avait calculé tous les points à partir de 0 :

# cette fois, nb_total1 et nb_cible1 sont remis à zéro
nb_cible2_verif = in_circle_plus(seed, 0, nb_total2, 0) 
print('nb_cible2_verif =', nb_cible2_verif)

nb_cible2_verif = 62833685

On obtient bien le même résultat. Mais le traitement est bien plus rapide avec advance(). Dans cet exemple, on a économisé pas moins de 7 huitièmes du temps de calcul.

# avec advance
%timeit in_circle_plus(seed, nb_total1, nb_total2, nb_cible1)
# sans advance
%timeit in_circle_plus(seed, 0, nb_total2, 0)

236 ms ± 1.08 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.94 s ± 39.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Cette méthode n’a bien sûr d’intérêt que si l’on veut reproduire une séquence de nombres aléatoires avec un seed déterminé.

En toute logique, on remarque que Generator(PCG(seed).advance(0)) génère le 1er terme de la suite, exactement comme si l’on n’avait pas utilisé la méthode advance().

Les trois instructions suivantes sont équivalentes et donnent les mêmes résultats :

rng1 = Generator(PCG64(42).advance(0))
rng2 = Generator(PCG64(42))
rng3 = np.random.default_rng(42)
for rng in (rng1, rng2, rng3):
    print(rng.random())

0.7739560485559633
0.7739560485559633
0.7739560485559633

Toutefois, NumPy évolue constamment et il n’est pas certain que le PCG64 restera le générateur par défaut dans les versions futures :

The default is currently PCG64 but this may change in future versions. As a convenience NumPy provides the default_rng function to hide these details. [doc]