02 · Séance type — simulation Monte-Carlo¶
Intuition
Quand un calcul est trop difficile à mener exactement, on le tire au sort. On simule le phénomène un grand nombre de fois, on observe la fréquence des résultats, et la loi des grands nombres garantit qu'elle converge vers la probabilité cherchée.
C'est brutal, c'est lent — l'erreur décroît seulement en \(\frac{1}{\sqrt n}\) — et c'est souvent la seule méthode disponible.
1. Le principe¶
Le schéma général
Pour estimer une quantité \(\theta = \mathbb{E}[g(X)]\) :
- simuler \(n\) réalisations indépendantes \(X_1,\dots,X_n\) ;
- calculer \(\hat\theta_n = \frac1n\sum_{i=1}^{n}g(X_i)\) ;
- estimer l'incertitude par \(\frac{\hat\sigma}{\sqrt n}\).
L'estimateur est sans biais (\(\mathbb{E}[\hat\theta_n]=\theta\)) et converge par la loi des grands nombres.
L'incertitude fait partie du résultat
Un résultat Monte-Carlo sans barre d'erreur est inexploitable. Par le TCL, l'intervalle de confiance à 95 % est
où \(\hat\sigma\) est l'écart-type empirique des \(g(X_i)\).
Coût de la précision : gagner un chiffre décimal exige 100 fois plus de tirages. C'est la limite fondamentale de la méthode.
2. Séance type — estimer \(\pi\)¶
Problème. Estimer \(\pi\) par tirage aléatoire.
Modélisation. On tire un point uniformément dans le carré \([0;1]^2\). La probabilité qu'il tombe dans le quart de disque unité est le rapport des aires :
Donc \(\pi = 4p\), et il suffit d'estimer \(p\) par la fréquence observée.
Hypothèses. Le générateur pseudo-aléatoire produit des couples \((U_1,U_2)\) uniformes et indépendants. C'est l'hypothèse critique — voir la critique en fin de séance.
import math
import random
def estimer_pi(n, rng):
dedans = 0
for _ in range(n):
x, y = rng.random(), rng.random()
if x * x + y * y <= 1.0:
dedans += 1
p = dedans / n
# ecart-type d'une Bernoulli, propage a pi = 4p
sigma_pi = 4 * math.sqrt(p * (1 - p) / n)
return 4 * p, 1.96 * sigma_pi
for n in (10**3, 10**5, 10**7):
rng = random.Random(42)
est, marge = estimer_pi(n, rng)
print(f"n={n:>9} pi ~ {est:.5f} +/- {marge:.5f} "
f"(erreur reelle {abs(est - math.pi):.5f})")
Résultats.
n= 1000 pi ~ 3.12800 +/- 0.10236 (erreur reelle 0.01359)
n= 100000 pi ~ 3.13728 +/- 0.01020 (erreur reelle 0.00431)
n= 10000000 pi ~ 3.14231 +/- 0.00102 (erreur reelle 0.00072)
Lecture des résultats
- L'erreur réelle est à chaque fois inférieure à la marge annoncée ✓ C'est le comportement attendu dans 95 % des cas.
- La marge est divisée par 10 quand \(n\) est multiplié par 100 — confirmation de la vitesse en \(\frac{1}{\sqrt n}\).
- Pour obtenir 6 décimales de \(\pi\), il faudrait \(n\approx10^{13}\) tirages. Monte-Carlo est une très mauvaise méthode pour calculer \(\pi\).
La critique qui compte
Cette séance a une valeur pédagogique, pas pratique. On connaît \(\pi\) à \(10^{-50}\) par des séries qui convergent en quelques dizaines de termes.
L'intérêt de Monte-Carlo n'apparaît que là où aucune autre méthode n'existe : intégrales en grande dimension, systèmes complexes, files d'attente non markoviennes.
3. Séance type — intégration en grande dimension¶
C'est le cas où Monte-Carlo est irremplaçable.
Problème. Calculer \(\displaystyle I = \int_{[0;1]^d}\cos\!\left(\sum_{i=1}^d x_i\right)\mathrm{d}\mathbf{x}\).
Solution exacte. Elle existe :
C'est précieux : elle sert de référence pour valider le code.
import cmath
import math
import random
def exact(d):
return ((cmath.exp(1j) - 1) / 1j) ** d
def monte_carlo(d, n, rng):
s = s2 = 0.0
for _ in range(n):
v = math.cos(sum(rng.random() for _ in range(d)))
s += v
s2 += v * v
moy = s / n
var = max(s2 / n - moy * moy, 0.0)
return moy, 1.96 * math.sqrt(var / n)
for d in (1, 5, 20):
rng = random.Random(3)
est, marge = monte_carlo(d, 200_000, rng)
ref = exact(d).real
print(f"d={d:2d} MC = {est:+.6f} +/- {marge:.6f} exact = {ref:+.6f}")
Résultats.
d= 1 MC = +0.841318 +/- 0.000607 exact = +0.841471
d= 5 MC = -0.649522 +/- 0.001630 exact = -0.649331
d=20 MC = -0.363023 +/- 0.002686 exact = -0.362095
Le point essentiel
La marge d'erreur ne dépend quasiment pas de \(d\). Elle passe de \(6\times10^{-4}\) en dimension 1 à \(2{,}7\times10^{-3}\) en dimension 20 : un facteur 4 seulement, alors que la dimension a été multipliée par 20.
Une méthode de quadrature classique, avec seulement 10 points par axe, demanderait \(10^{20}\) évaluations en dimension 20 — irréalisable. Monte-Carlo en fait \(2\times10^5\).
C'est la seule méthode connue dont le coût est indépendant de la dimension. Voir l'exercice 9 du chapitre Probabilités 10.
4. Séance type — simulation d'un système¶
Problème. Un guichet, arrivées poissonniennes d'intensité \(\lambda\), temps de service exponentiels de taux \(\mu\). Quel est le temps d'attente moyen ?
Pourquoi simuler. La formule M/M/1 existe (\(W = \frac{\rho}{\mu(1-\rho)}\)), mais elle cesse d'exister dès qu'on change une hypothèse : deux guichets avec des vitesses différentes, clients impatients, priorités. La simulation, elle, s'adapte.
import math
import random
def simuler_file(lam, mu, n_clients, rng):
"""File M/M/1. Renvoie le temps d'attente moyen."""
t_arrivee = 0.0
t_libre = 0.0 # instant ou le guichet se libere
attentes = []
for _ in range(n_clients):
t_arrivee += -math.log(rng.random()) / lam # inter-arrivee ~ E(lambda)
debut = max(t_arrivee, t_libre)
attentes.append(debut - t_arrivee)
service = -math.log(rng.random()) / mu # service ~ E(mu)
t_libre = debut + service
return sum(attentes) / len(attentes)
lam, mu = 20.0, 25.0 # clients par heure
rho = lam / mu
theorique = rho / (mu * (1 - rho))
rng = random.Random(7)
simule = simuler_file(lam, mu, 500_000, rng)
print(f"rho = {rho}")
print(f"W theorique = {theorique*60:.2f} min")
print(f"W simule = {simule*60:.2f} min")
Résultats.
rho = 0.8
W theorique = 9.60 min
W simule = 9.79 min
La simulation valide la formule — et réciproquement
L'écart de 2 % est cohérent avec l'incertitude de simulation. La concordance valide simultanément le code et la compréhension du modèle.
C'est le geste méthodologique central : tester le simulateur sur un cas où la réponse est connue, avant de l'utiliser sur un cas où elle ne l'est pas.
Extension immédiate : en changeant deux lignes — service à durée constante au lieu d'exponentielle —, on obtient une file M/D/1, dont le temps d'attente est théoriquement la moitié. C'est vérifiable en trente secondes, et cela donne une conclusion opérationnelle : réduire la variabilité du service vaut autant qu'accélérer le service.
5. Réduction de variance¶
Quand la précision plafonne, on peut réduire \(\sigma\) plutôt qu'augmenter \(n\).
Trois techniques
| Technique | Principe | Gain typique |
|---|---|---|
| Variables antithétiques | Utiliser \(U\) et \(1-U\) par paire | ×2 à ×10 |
| Échantillonnage stratifié | Découper le domaine, tirer dans chaque strate | ×2 à ×100 |
| Échantillonnage préférentiel | Tirer davantage là où \(g\) est grande | ×10 à ×10⁶ |
Variables antithétiques — la plus simple. Si \(g\) est monotone, \(g(U)\) et \(g(1-U)\) sont négativement corrélées, donc
Le gain vient directement de la formule de la variance d'une somme, chapitre Probabilités 10.
Exercices¶
★★ Exercice 1. Estimer \(\displaystyle\int_0^1\sqrt{1-x^2}\,\mathrm{d}x\) par Monte-Carlo, avec intervalle de confiance. Comparer à la valeur exacte \(\frac\pi4\).
★★ Exercice 2. Estimer par simulation la probabilité qu'en lançant 5 dés on obtienne au moins une paire. Comparer au calcul exact.
★★★ Exercice 3. Implémenter les variables antithétiques pour l'intégrale de l'exercice 1 et mesurer le gain de variance.
★★★ Exercice 4. Simuler le problème du collectionneur de coupons pour \(n=50\) et comparer à la valeur théorique \(nH_n \approx 225\) (chapitre Probabilités 07, exercice 9).
★★★★ Exercice 5. Simuler une file M/D/1 (service de durée constante \(\frac1\mu\)) et vérifier que le temps d'attente vaut la moitié de celui de la file M/M/1 de mêmes paramètres.
Rédiger la critique : pourquoi ce résultat, et qu'implique-t-il en pratique ?
★★★★ Exercice 6. Estimer par Monte-Carlo la probabilité qu'une matrice \(3\times3\) à coefficients uniformes sur \([-1;1]\) ait un déterminant positif.
a) Prédire le résultat par un argument de symétrie avant de simuler. b) Simuler et vérifier. c) Que se passe-t-il pour une matrice \(n\times n\) ?
Corrigés¶
Corrigé — Exercice 1
import math
import random
def integrale_mc(g, n, rng):
s = s2 = 0.0
for _ in range(n):
v = g(rng.random())
s += v
s2 += v * v
moy = s / n
var = max(s2 / n - moy * moy, 0.0)
return moy, 1.96 * math.sqrt(var / n)
g = lambda x: math.sqrt(1 - x * x)
for n in (10**4, 10**6):
rng = random.Random(1)
est, marge = integrale_mc(g, n, rng)
print(f"n={n:>8} I = {est:.5f} +/- {marge:.5f} exact = {math.pi/4:.5f}")
Résultats.
n= 10000 I = 0.78555 +/- 0.00440 exact = 0.78540
n= 1000000 I = 0.78525 +/- 0.00044 exact = 0.78540
L'erreur réelle (\(1{,}5\times10^{-4}\) puis \(1{,}5\times10^{-4}\)) reste dans la marge annoncée ✓ et la marge est divisée par 10 quand \(n\) est multiplié par 100 ✓
Corrigé — Exercice 2
Calcul exact d'abord. « Au moins une paire » est le complémentaire de « 5 valeurs toutes distinctes » :
import random
def au_moins_une_paire(rng, n):
succes = 0
for _ in range(n):
des = [rng.randint(1, 6) for _ in range(5)]
if len(set(des)) < 5:
succes += 1
return succes / n
rng = random.Random(5)
p = au_moins_une_paire(rng, 200_000)
exact = 1 - (6 * 5 * 4 * 3 * 2) / 6**5
print(f"simule = {p:.5f} exact = {exact:.5f}")
Résultat : simule = 0.90583 exact = 0.90741
Écart de \(1{,}6\times10^{-3}\) — compatible avec la marge \(1{,}96\sqrt{\frac{0{,}906\times0{,}094}{200000}} = 1{,}3\times10^{-3}\), à la limite de l'intervalle : c'est l'un des 5 % de tirages qui en sortent.
Corrigé — Exercice 3
import math
import random
g = lambda x: math.sqrt(1 - x * x)
def mc_simple(n, rng):
vals = [g(rng.random()) for _ in range(n)]
m = sum(vals) / n
var = sum((v - m) ** 2 for v in vals) / n
return m, var
def mc_antithetique(n, rng):
# n/2 paires (U, 1-U) : meme nombre d'appels a g
vals = []
for _ in range(n // 2):
u = rng.random()
vals.append((g(u) + g(1 - u)) / 2)
m = sum(vals) / len(vals)
var = sum((v - m) ** 2 for v in vals) / len(vals)
return m, var
n = 200_000
rng = random.Random(11)
m1, v1 = mc_simple(n, rng)
rng = random.Random(11)
m2, v2 = mc_antithetique(n, rng)
# variance de l'estimateur final (et non de la variable)
var_est_simple = v1 / n
var_est_anti = v2 / (n // 2)
print(f"simple : I = {m1:.6f} var(estimateur) = {var_est_simple:.3e}")
print(f"antithetique : I = {m2:.6f} var(estimateur) = {var_est_anti:.3e}")
print(f"gain de variance : x{var_est_simple/var_est_anti:.2f}")
Résultat typique :
simple : I = 0.785484 var(estimateur) = 2.486e-07
antithetique : I = 0.785272 var(estimateur) = 6.886e-08
gain de variance : x3.61
La variance est divisée par environ 3,6, à nombre d'appels à \(g\) égal. C'est équivalent à multiplier \(n\) par 3,6 — gratuitement.
Pourquoi ça marche : \(g\) est décroissante sur \([0;1]\), donc \(g(U)\) et \(g(1-U)\) sont fortement anticorrélées (\(\rho\approx-0{,}5\)), et la formule de la variance d'une somme donne le gain.
Corrigé — Exercice 4
import math
import random
def collectionneur(n, rng):
vus = set()
achats = 0
while len(vus) < n:
vus.add(rng.randrange(n))
achats += 1
return achats
n = 50
rng = random.Random(2026)
essais = [collectionneur(n, rng) for _ in range(20_000)]
moyenne = sum(essais) / len(essais)
theorique = n * sum(1 / k for k in range(1, n + 1))
ecart_type = math.sqrt(sum((x - moyenne)**2 for x in essais) / len(essais))
print(f"simule = {moyenne:.2f} (ecart-type {ecart_type:.1f})")
print(f"theorique = {theorique:.2f}")
print(f"mediane = {sorted(essais)[len(essais)//2]}")
Résultat typique :
simule = 224.91 (ecart-type 61.9)
theorique = 224.96
mediane = 214
Concordance à 0,1 % ✓
Observation importante : l'écart-type est de 61 achats pour une moyenne de 225 — soit 27 %. La distribution est très dispersée et asymétrique (médiane 214 < moyenne 225).
Annoncer « il faut 225 achats » est donc trompeur : il en faut entre 150 et 350 selon la chance. C'est exactement le genre de nuance que la critique doit apporter.
Corrigé — Exercice 5
import math
import random
def simuler(lam, mu, n_clients, rng, service_constant=False):
t_arrivee = t_libre = 0.0
attentes = []
for _ in range(n_clients):
t_arrivee += -math.log(rng.random()) / lam
debut = max(t_arrivee, t_libre)
attentes.append(debut - t_arrivee)
if service_constant:
s = 1.0 / mu
else:
s = -math.log(rng.random()) / mu
t_libre = debut + s
return sum(attentes) / len(attentes)
lam, mu = 20.0, 25.0
rho = lam / mu
rng = random.Random(7)
mm1 = simuler(lam, mu, 500_000, rng)
rng = random.Random(7)
md1 = simuler(lam, mu, 500_000, rng, service_constant=True)
print(f"M/M/1 simule = {mm1*60:.2f} min (theorie {rho/(mu*(1-rho))*60:.2f})")
print(f"M/D/1 simule = {md1*60:.2f} min (theorie {rho/(2*mu*(1-rho))*60:.2f})")
print(f"rapport = {mm1/md1:.3f}")
Résultat typique :
M/M/1 simule = 9.79 min (theorie 9.60)
M/D/1 simule = 4.75 min (theorie 4.80)
rapport = 2.062
Le rapport vaut 2,06, soit 2 à 3 % près ✓
Critique — pourquoi ce résultat. La formule de Pollaczek-Khintchine donne, pour une file M/G/1 avec un service de moyenne \(\frac1\mu\) et de coefficient de variation \(c_s\) :
Pour un service exponentiel, \(c_s = 1\) et le facteur vaut 1. Pour un service déterministe, \(c_s = 0\) et le facteur vaut \(\frac12\).
L'attente est donc directement proportionnelle à \(1+c_s^2\) : c'est la variabilité du service, pas seulement sa vitesse, qui crée l'attente.
Implication pratique. Réduire la variance du temps de traitement — en normalisant les procédures, en écartant les cas atypiques vers une file dédiée — divise l'attente par deux sans accélérer quoi que ce soit.
C'est un résultat contre-intuitif et opérationnel : dans un centre d'appels, séparer les demandes simples des demandes complexes améliore l'attente moyenne des deux catégories.
Corrigé — Exercice 6
a) Prédiction par symétrie.
Échanger deux lignes d'une matrice change le signe du déterminant (chapitre Algèbre 07).
Or, si \(\mathbf{A}\) a ses coefficients tirés i.i.d. selon une loi symétrique (ici uniforme sur \([-1;1]\)), la matrice obtenue en échangeant deux lignes a exactement la même loi.
L'application \(\mathbf{A}\mapsto\mathbf{A}'\) (échange des lignes 1 et 2) est donc une bijection préservant la loi et changeant le signe du déterminant. D'où
Et \(P(\det=0)=0\) (l'ensemble des matrices singulières est de mesure nulle).
b) Vérification.
import random
def det3(m):
a, b, c = m[0]
d, e, f = m[1]
g, h, i = m[2]
return a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)
rng = random.Random(99)
n = 200_000
positifs = 0
for _ in range(n):
m = [[rng.uniform(-1, 1) for _ in range(3)] for _ in range(3)]
if det3(m) > 0:
positifs += 1
p = positifs / n
marge = 1.96 * (p * (1 - p) / n) ** 0.5
print(f"P(det>0) = {p:.5f} +/- {marge:.5f}")
Résultat : P(det>0) = 0.49838 +/- 0.00219
L'intervalle contient 0,5 ✓ La prédiction théorique est confirmée.
c) Cas \(n\times n\). L'argument de symétrie ne dépend pas de \(n\) — il suffit que \(n\geqslant2\) pour qu'un échange de lignes soit possible.
Pour \(n=1\), en revanche, \(\det\mathbf{A} = a_{11}\) et \(P(\det>0)=\frac12\) également — par symétrie de la loi uniforme cette fois.
La leçon méthodologique
Prédire avant de simuler. Une simulation qui confirme une prédiction valide les deux ; une simulation sans prédiction ne valide rien — on ne saurait pas distinguer un vrai résultat d'un bug.
C'est la règle la plus importante de tout ce chapitre.
Chapitre suivant : Séance type — ajustement de données.