En prépa ECG, simuler une loi de probabilité classique ne pose pas de difficulté : Python dispose de commandes toutes faites pour la loi normale, la loi exponentielle ou la loi binomiale. Mais que faire lorsque la loi que l’on cherche à simuler est définie par une densité quelconque, sans commande dédiée ? C’est précisément ce que résout la méthode de rejet, aussi appelée méthode d’acceptation-rejet. Elle permet de simuler n’importe quelle variable aléatoire à densité à partir de simples tirages uniformes, en s’appuyant sur un raisonnement probabiliste élégant entièrement au programme. Dans cet article, nous allons présenter le principe de la méthode et le justifier rigoureusement, étudier le choix optimal de la densité instrumentale et son impact sur l’efficacité, puis implémenter la méthode en Python sur deux exemples concrets.
Principe et justification de la méthode
Le problème posé
On se donne une densité de probabilité \(f\) sur un intervalle \([a, b]\), que l’on souhaite simuler. On suppose que \(f\) est bornée sur \([a, b]\), c’est-à-dire qu’il existe \(M > 0\) tel que \(f(x) \leq M\) pour tout \(x \in [a, b]\). La méthode de rejet repose sur l’observation géométrique suivante : simuler une variable aléatoire de densité \(f\) revient à tirer un point uniformément sous la courbe de \(f\) et à ne conserver que son abscisse.
Plus précisément, on considère le rectangle \([a, b] \times [0, M]\) qui contient entièrement le graphe de \(f\). Si l’on tire un point \((X, U)\) uniformément dans ce rectangle, la probabilité que ce point tombe sous la courbe de \(f\) est :
\[\frac{\int_a^b f(x)\,dx}{M(b-a)} = \frac{1}{M(b-a)}\]
puisque \(f\) est une densité. Conditionnellement à l’événement \(\{U \leq f(X)\}\), la variable \(X\) suit exactement la loi de densité \(f\). C’est ce conditionnement qui constitue le cœur de la méthode.
L’algorithme
L’algorithme est le suivant. On tire \(X\) selon une loi uniforme sur \([a, b]\) et \(U\) selon une loi uniforme sur \([0, M]\), indépendamment. Si \(U \leq f(X)\), on accepte \(X\) comme réalisation de la loi cible. Sinon, on rejette et on recommence.
Voici un graph illustrant la méthode :
Tu pourras remarquer des similitudes avec la méthode de Monte-Carlo, puisqu’il y a aussi des simulations de lois uniformes en jeu.
Justification rigoureuse
Montrons que \(X\) conditionné à l’acceptation suit bien la loi de densité \(f\). Notons \(A\) l’événement d’acceptation \(\{U \leq f(X)\}\). Pour tout borélien \(B \subset [a, b]\) :
\[P(X \in B \mid A) = \frac{P(X \in B, U \leq f(X))}{P(A)}\]
Or, \(X\) et \(U\) étant indépendants et uniformément distribués sur \([a,b]\) et \([0, M]\) respectivement :
\[P(X \in B, U \leq f(X)) = \int_B \frac{1}{b-a} \cdot \frac{f(x)}{M}\,dx = \frac{1}{M(b-a)} \int_B f(x)\,dx\]
Et \(P(A) = \frac{1}{M(b-a)}\), comme vu dans la première partie. On obtient donc :
\[P(X \in B \mid A) = \frac{\frac{1}{M(b-a)} \int_B f(x)\,dx}{\frac{1}{M(b-a)}} = \int_B f(x)\,dx\]
Ce qui est bien la probabilité que la variable de densité \(f\) appartienne à \(B\). La justification est complète.
Efficacité et choix de la constante M
Le taux d’acceptation
La probabilité d’accepter un tirage est \(P(A) = \frac{1}{M(b-a)}\). Plus \(M\) est grand, plus cette probabilité est faible, et plus il faut de tirages pour obtenir une réalisation acceptée. En moyenne, le nombre de tirages nécessaires pour obtenir une réalisation est :
\[\frac{1}{P(A)} = M(b-a)\]
Il est donc crucial de choisir \(M\) aussi petit que possible, tout en garantissant \(f(x) \leq M\) pour tout \(x \in [a, b]\). Le choix optimal est \(M = \sup_{x \in [a,b]} f(x)\), c’est-à-dire le maximum de la densité cible sur son support.
Extension à un support non borné
Attention, cependant ! La méthode telle que présentée suppose \(f\) à support borné. Lorsque le support est \(\mathbb{R}\) entier (comme pour la loi normale), on ne peut pas tirer \(X\) uniformément sur \(\mathbb{R}\). La solution est de remplacer la loi uniforme par une densité instrumentale \(g\) à support identique à celui de \(f\) et de chercher une constante \(M\) telle que \(f(x) \leq M \cdot g(x)\) pour tout \(x\).
L’algorithme devient alors : tirer \(X\) selon \(g\) et \(U\) selon une loi uniforme sur \([0,1]\) et accepter si \(U \leq \frac{f(X)}{M \cdot g(X)}\). Le taux d’acceptation vaut alors \(\frac{1}{M}\), indépendamment du support. Plus \(g\) est proche de \(f\), plus \(M\) peut être choisi petit, et plus la méthode est efficace.
Voilà pour cette extension à support non borné que nous ne développerons pas davantage, puisque cette notion s’écarte du cœur de cet article.
Implémentation en Python
Exemple 1 : simuler la loi bêta
La loi bêta de paramètres \(\alpha = 2\) et \(\beta = 3\) a pour densité sur \([0, 1]\) :
\[f(x) = \frac{x^{\alpha – 1}(1-x)^{\beta – 1}}{B(\alpha, \beta)} = 12\,x(1-x)^2\]
Son maximum sur \([0,1]\) est atteint en \(x^* = \frac{\alpha – 1}{\alpha + \beta – 2} = \frac{1}{3}\), ce qui donne \(M = f\!\left(\frac{1}{3}\right) = 12 \cdot \frac{1}{3} \cdot \frac{4}{9} = \frac{16}{9} \approx 1{,}78\). Le taux d’acceptation théorique est \(\frac{1}{M \cdot (1-0)} = \frac{1}{M} \approx 56\%\).
On obtient alors le graph suivant : 
Analysons ce script
On définit d’abord la densité (f) de la loi bêta et la constante \(M = \frac{16}{9} \), son maximum sur ([0,1]). La boucle while tourne jusqu’à avoir accumulé (n) réalisations acceptées : à chaque itération, on tire (X) uniformément sur ([0,1]) et (U) uniformément sur ([0, M]) et on n’accepte (X) que si \(U \leq f(X)\), c’est-à-dire si le point ((X, U)) tombe sous la courbe de (f).
Le compteur nb_tirages permet de calculer empiriquement le taux d’acceptation en fin de simulation, que l’on compare à la valeur théorique \(\frac{1}{M} = \frac{9}{16} \approx 56\%\). Enfin, on superpose l’histogramme normalisé des réalisations acceptées à la densité théorique : si les deux coïncident visuellement, c’est la confirmation que l’algorithme simule bien la loi cible.
À l’évidence, ce long script permet de te montrer visuellement que cette méthode fonctionne bien. Mais en pratique, il suffirait de s’arrêter à la ligne 20 du script précédent pour simuler la loi bêta.
Exemple 2 : simuler la loi normale tronquée
On veut simuler la loi normale \(\mathcal{N}(0,1)\) tronquée sur \([-3, 3]\), dont la densité est :
\[f(x) = \frac{1}{\sqrt{2\pi}}\,e^{-x^2/2} \cdot \frac{1}{P(-3 \leq Z \leq 3)}\]
On prend comme densité instrumentale la loi uniforme sur \([-3, 3]\), soit \(g(x) = \frac{1}{6}\). La constante \(M\) est le rapport maximal \(\frac{f(x)}{g(x)}\), atteint en \(x = 0\) :
\[M = \frac{f(0)}{g(0)} = \frac{6}{\sqrt{2\pi}} \approx 2{,}39\]
Le taux d’acceptation théorique est \(\frac{1}{M} \approx 42\%\), ce que le script ci-dessous confirme empiriquement.
return np.array(resultats) renverra alors les réalisations de la loi normale centrée réduite tronquée sur \([-3;3]\).
Ces deux exemples illustrent bien la démarche générale : identifier le maximum de \(f\) pour fixer \(M\), puis appliquer le test d’acceptation. La seule vraie difficulté pratique est le calcul de \(M\), qui peut nécessiter une étude de variations préalable.
Conclusion
En définitive, la méthode de rejet est un algorithme universel pour simuler n’importe quelle loi à densité, à condition de pouvoir borner cette densité. Elle s’appuie entièrement sur des outils au programme : loi uniforme, conditionnement, calcul intégral pour la justification et étude de variations pour optimiser \(M\).
Bien que hors programme ECG au sens strict, elle peut très bien apparaître dans un sujet de Parisiennes comme extension naturelle de la simulation de lois classiques. Savoir la coder rapidement et expliquer le raisonnement derrière le test d’acceptation peut rapporter des points décisifs.
Il n’y a plus qu’à croiser les doigts pour qu’un tel thème tombe aux concours !
Tu peux retrouver ici le méga-répertoire qui contient toutes les annales de concours et les corrigés. Tu peux également accéder ici à toutes nos autres ressources mathématiques !
