Perspectives historiques
{{
}}
Nous allons présenter ici quelques éléments historiques sur les méthodes de Monte-Carlo, dont les prémisses remontent au XVIIIème siècle.
L’aiguille de Buffon
Georges-Louis Leclerc, Comte de Buffon1 proposa en 1733 une méthode qui s’avéra être utile pour estimer la valeur de \pi. On désigne de nos jours cette expérience sous le nom de l’aiguille de Buffon. C’est l’une des premières méthodes de Monte-Carlo référencée dans la littérature (la source du texte est disponible ici sur le site de la BNF).
1 Georges-Louis Leclerc, Comte de Buffon: (1707-1788) naturaliste, mathématicien et industriel français du siècle des Lumières 
La question initiale (simplifiée ici) posée par Buffon était la suivante: une aiguille de taille 1 tombe sur un parquet composé de lattes de largeur 1: quelle est alors la probabilité P que l’aiguille croise une ligne de la trame du parquet ?
Le contexte original était dans celui d’un jeu à deux joueurs: un joueur parie sur le fait que l’aiguille croise une ligne de la trame du parquet, l’autre sur le fait que l’aiguille ne croise pas une ligne de la trame du parquet. L’enjeu est alors de calculer la probabilité de succès de chacun des joueurs, et de voir si le jeu est équilibré ou non.
Voilà brièvement la question que s’est posée Buffon en 1733. La réponse est donnée par la formule suivante, qui montre que le jeu qu’il propose n’est pas équilibré:
P = \frac{2}{\pi} \approx 0.6366 \enspace. Une preuve de ce résultat sera donnée ci-dessous.
L’idée sous-jacente de Buffon est que si l’on répète cette expérience un grand nombre de fois, on peut approché la quantité P numériquement, par exemple en proposant un estimateur \hat{P}_n qui compte la proportion de chevauchement après avoir fait n répétition des lancers. Pour estimer \pi, il ne restera donc plus qu’à évaluer \frac{2}{\hat{P}_n}.
On peut faire cette expérience dans le monde réel (c’est un peu long pour n grand!), mais on peut aussi utiliser une méthode numérique pour cela. Il s’agit alors de tirer aléatoirement la position du centre de l’aiguille, puis de tirer aussi de manière aléatoire son angle de chute. On teste à la fin si l’aiguille croise une ligne de la trame du parquet ou non, et on recommence l’expérience un grand nombre de fois.
Cette méthode est donnée ci-dessous, avec un exemple interactif généré en Python.
Code
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
rng = np.random.default_rng(44)
n_samples = 200
xmax = 14.499999
xmin = -xmax
# Create the needles
centers_x = rng.uniform(xmin, xmax, n_samples)
angles = rng.uniform(0, 2 * np.pi, n_samples)
centers_y = rng.uniform(-2, 2, n_samples)
# Compute the right borders of the needles
borders_right = np.zeros((n_samples, 2))
borders_right[:, 0] = centers_x + np.cos(angles) / 2
borders_right[:, 1] = centers_y + np.sin(angles) / 2
# Compute the left borders of the needles
borders_left = np.zeros((n_samples, 2))
borders_left[:, 0] = centers_x + np.cos(angles + np.pi) / 2
borders_left[:, 1] = centers_y + np.sin(angles + np.pi) / 2
centers_x_round = np.round(centers_x)
overlap = (borders_right[:, 0] - centers_x_round) * (
borders_left[:, 0] - centers_x_round
) < 0
overlap = np.where(overlap, 1, 0)
n_overlap = int(np.sum(overlap))
# Check if the needles cross a line
borders_red = np.empty((3 * n_overlap, 2), dtype=object)
borders_red.fill(None)
borders_red[::3, :] = borders_right[overlap == 1]
borders_red[1::3, :] = borders_left[overlap == 1]
borders_blue = np.empty((3 * (n_samples - n_overlap), 2), dtype=object)
borders_blue.fill(None)
borders_blue[::3, :] = borders_right[overlap == 0]
borders_blue[1::3, :] = borders_left[overlap == 0]
overlaps = np.empty((3 * n_samples), dtype=object)
overlaps.fill(None)
overlaps[::3] = overlap
overlaps[1::3] = overlap
overlaps[2::3] = overlap
idx_red = np.cumsum(overlaps)
idx_blue = np.cumsum(1 - overlaps)
# Create subplots with 2 rows and 1 column with ratio x / y of 10
fig = make_subplots(rows=2, cols=1, vertical_spacing=0.1, row_heights=[2, 1])
# Use a loop to plot vertical lines equation "y=c" for integer values c in [-2, -1, 0, 1, 2]
for i in range(int(np.round(xmin)), int(np.round(xmax)) + 1):
fig.add_shape(
type="line",
y0=-3,
x0=i,
y1=3,
x1=i,
line=dict(
color="black",
width=2,
),
row=1,
col=1,
)
color = np.where(overlaps, 1.0, 0.0)
n_samples_array = np.arange(1, n_samples + 1)
pi_estimate = 2 / (np.cumsum(overlap) / n_samples_array)
t = n_samples
fig.update_layout(
template="simple_white",
xaxis=dict(range=[xmin, xmax], constrain="domain", showgrid=False),
yaxis_scaleanchor="x",
xaxis_visible=False,
yaxis_visible=False,
)
for i in range(3, t):
fig.add_trace(
go.Scatter(
x=borders_red[: idx_red[3 * i] + 1, 0],
y=borders_red[: idx_red[3 * i] + 1, 1],
mode="lines",
line=dict(width=2),
marker=dict(color="red"),
name="Avec intersection",
visible=False,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=borders_blue[: idx_blue[3 * i] + 1, 0],
y=borders_blue[: idx_blue[3 * i] + 1, 1],
mode="lines",
line=dict(width=2),
marker=dict(color="darkblue"),
name="Sans intersection",
visible=False,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=n_samples_array[:i],
y=pi_estimate[:i],
mode="lines",
line=dict(width=1),
marker=dict(color="red"),
showlegend=False,
visible=False,
),
row=2,
col=1,
)
fig.add_annotation(
dict(
x=1.01,
y=0.14,
xref="paper",
yref="paper",
text="Estimation de pi",
showarrow=False,
font=dict(color="red"),
)
)
fig.add_annotation(
dict(x=-0.04, y=0.19, xref="paper", yref="paper", text="pi", showarrow=False)
)
fig.update_xaxes(title_text="Nombre d'aiguilles tirées", row=2, col=1)
fig.update_layout(
template="none",
xaxis2=dict(showgrid=True, zeroline=True, zerolinewidth=1, range=[0, n_samples]),
yaxis2=dict(showgrid=True, zeroline=True, zerolinewidth=1, range=[0, 6]),
)
# plot a dash line at y=pi
fig.add_shape(
type="line",
y0=np.pi,
x0=0,
y1=np.pi,
x1=n_samples,
line=dict(
color="black",
width=1,
dash="dashdot",
),
row=2,
col=1,
)
fig.data[10 * 3].visible = True
fig.data[10 * 3 + 1].visible = True
fig.data[10 * 3 + 2].visible = True
steps = []
for i in range(len(fig.data) // 3):
step = dict(
label=str(i + 4),
method="update",
args=[
{"visible": [False] * len(fig.data)},
{
"title": "Estimation avec "
+ str(i + 4)
+ f" aiguilles: pi = {pi_estimate[i]:.4f}"
},
],
)
step["args"][0]["visible"][3 * i] = True
step["args"][0]["visible"][3 * i + 1] = True
step["args"][0]["visible"][3 * i + 2] = True
steps.append(step)
slider = dict(
active=0,
currentvalue={"prefix": "Nombre d'aiguilles: "},
pad={"t": 50},
y=-0.32,
steps=steps,
)
fig.update_layout(legend=dict(x=0.5, y=0.31, xanchor='center', yanchor='bottom'))
fig.update_layout(sliders=[slider])
fig.show()On va fournir ici le calcul de la probabilité P. Pour cela on aura besoin de quelques éléments décrits dans le dessin ci-dessous.
- x : distance entre le centre de l’aiguille et la ligne de la trame du parquet la plus proche
- \theta : angle entre l’aiguille et la ligne de la trame du parquet la plus proche
- 1 : longueur de l’aiguille (et donc la demi longueur est \frac{1}{2})
- \frac{1}{2}\sin(\theta) : distance entre l’extrémité de l’aiguille et la ligne de la trame du parquet la plus proche
Avec les éléments ci-dessus, on voit qu’il y a chevauchement si et seulement si: \frac{1}{2}\sin(\theta) \geq x.
Maintenant par des arguments de symétrie on voit qu’on peut se restreindre à \theta \in [0, \frac{\pi}{2}], et à x \in [0, \frac{1}{2}]. Les lois de générations des variables aléatoires X et \Theta sont les suivantes:
- X \sim \mathcal{U}([0, \frac{1}{2}]), de densité f_X(x) = 2 {1\hspace{-3.8pt} 1}_{[0, \frac{1}{2}]}(x)
- \Theta \sim \mathcal{U}([0, \frac{\pi}{2}]) de densité f_\Theta(\theta) = \frac{2}{\pi} {1\hspace{-3.8pt} 1}_{[0, \frac{\pi}{2}]}(\theta)
De plus on suppose que X et \Theta sont indépendantes.
Maintenant pour calculer la probabilité P on procède comme suit:
\begin{align*} P & = \mathbb{P}\left(\frac{1}{2}\sin(\Theta) \geq X\right) \\ & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{\frac{1}{2}\sin(\theta) \geq x\}} f_{\Theta}(\theta) f_X(x) d\theta dx \quad (\text{par indépendance})\\ & = \int_{\mathbb{R}} \int_{\mathbb{R}} {1\hspace{-3.8pt} 1}_{\{\frac{1}{2}\sin(\theta) \geq x\}} \frac{2}{\pi} {1\hspace{-3.8pt} 1}_{[0, \frac{\pi}{2}]}(\theta) \cdot 2 {1\hspace{-3.8pt} 1}_{[0, \frac{1}{2}]}(x) d\theta dx \\ & = \int_{0}^{\frac{\pi}{2}} \int_{0}^{\frac{1}{2}\sin(\theta)} \frac{4}{\pi} dx d\theta \\ & = \frac{4}{\pi} \int_{0}^{\frac{\pi}{2}} {\frac{1}{2}\sin(\theta)} d\theta\\ & = \frac{2}{\pi} \Big[ -\cos(\theta)\Big]_{0}^{\frac{\pi}{2}} \\ & = \frac{2}{\pi} \enspace. \end{align*}
Les spaghettis de Buffon
Une démonstration alternative est donnée dans le livre de Borel et Deltheil (1923), inspirée de celle Barbier (1860), et utilise concept de spaghettis de Buffon (🇬🇧 Buffon’s noodle, cf. https://en.wikipedia.org/wiki/Buffon%27s_noodle). Elle utilise l’approche suivante. Au lieu de calculer la probabilité d’intersection entre une aiguille et les rainures d’un parquet, nous allons déterminer le nombre moyen d’intersections entre une courbe et les rainures du parquet. Nous démontrerons que ce nombre moyen est directement proportionnel à la longueur de la courbe et inversement proportionnel à la distance moyenne entre les rainures du parquet (nous prendrons cette distance égale à 1 pour simplifier).
Pour ce faire, nous allons d’abord considérer des lignes polygonales, pour lesquelles l’espérance d’intersection se calcule en sommant les espérances sur chaque segment. Nous étendrons ensuite ce résultat par continuité aux courbes régulières. La constante de proportionnalité de cette relation sera notée \alpha.
Nous allons maintenant montrer que \alpha = \frac{2}{\pi}. Pour cela, prenons un cercle de rayon 2, donc de longueur 4\pi. Si ce cercle est tangent aux lattes du parquet, il y a deux points d’intersection. Sinon, il y a également deux points d’intersection. Ainsi, le nombre moyen d’intersection est de 2, ce qui nous donne l’équation 2 = \alpha \times 4\pi. En résolvant cette équation, nous trouvons \alpha = \frac{2}{\pi}.
Ce résultat nous donne la probabilité d’intersection entre une aiguille et les lattes du parquet lorsque la taille de l’aiguille est inférieure à l’écartement des lattes.
Méthode de Monte-Carlo
La méthode de Monte-Carlo, est une méthode de calcul numérique qui consiste à utiliser des nombres aléatoires pour résoudre des problèmes déterministes. Elle est utilisée dans de nombreux domaines, comme la physique, la chimie, la biologie, la finance, ou encore l’apprentissage automatique. Cette méthode basée sur la loi des grands nombres a été mis au point à Los Alamos, dans le cadre du projet Manhattan (dont l’objectif était le développement du nucléaire civil et militaire) par un groupe de scientifiques dont les plus connus sont: John von Neumann2, Nicholas Metropolis3 ou encore Stanisław Ulam4
2 John von Neumann: (1903-1957) mathématicien et physicien américano-hongrois, un des pères de l’informatique. 
3 Nicholas Metropolis: (1915-1999), physicien gréco-américain, est des initiateurs de la méthode de Monte Carlo et du recuit simulé 
4 Stanisław Ulam: (1909-1984) 
Dans le cadre du projet Manhattan, il s’agissait de calculer des intégrales de manière numérique pour modéliser l’évolution de particules, en utilisant des nombres aléatoires.
Eckhardt (1987) donne un bref aperçu historique, et mentionne les premières description de la méthode du rejet et de la méthode de l’inversion dans des lettres entre Von Neumann et Ulam datant de 1947. Ulam aurait une l’idée d’utiliser de telles méthodes pour résoudre le jeu du solitaire lors d’un séjour à l’hôpital en 1946, et éviter ainsi de faire des calculs combinatoires fastidieux. Rapidement, la possibilité d’appliquer cette approche pour des calculs en physique mathématiques (diffusion des neutrons notamment) lui serait apparue prometteuse. Le développement de l’informatique naissante allait permettre une mise en oeuvre pratique de ces idées, et c’est ainsi que la méthode de Monte-Carlo est née. Le nom Monte-Carlo est lui venu du besoin de confidentialité du projet, et provient du nom de la ville de Monte-Carlo, connue pour ses jeux de hasard, où l’oncle de Stanisław Ulam aimait se rendre pour assouvir sa soif de jeu. Ce serait N. Metropolis qui aurait proposé ce nom (cf. Metropolis 1987):
It was at that time that I suggested an obvious name for the statistical method—a suggestion not unrelated to the fact that Stan had an uncle who would borrow money from relatives because he “just had to go to Monte Carlo”.
Autres méthodes stochastiques populaires
Méthode d’Hasting-Metropolis
L’algorithme de Hasting-Metropolis est une méthode MCMC (🇬🇧: Monte Carlo Markov Chains) dont le but est d’obtenir un échantillonnage aléatoire d’une distribution de probabilité quand l’échantillonnage direct en est difficile (en particulier en grande dimension)
Un avantage est qu’il ne requiert la connaissance de loi de densité qu’à constante multiplicative près.
Recuit simulé
Le recuit simulé ( 🇬🇧 simulated annealing) est une méthode (empirique) d’optimisation, inspirée d’un processus, le recuit, utilisé en métallurgie. On alterne dans cette dernière des cycles de refroidissement lent et de réchauffage (recuit) qui ont pour effet de minimiser l’énergie du matériau. Cette méthode est transposée en optimisation pour trouver les extrema d’une fonction.