Simulation

Dans ce chapitre on se demande comment simuler en pratique des variables aléatoires i.i.d. L’idée est de commencer par le cas de variables aléatoires de loi uniforme et d’en déduire les autres lois.

Variables aléatoires uniformes

On rappelle qu’une variable aléatoire U suit une loi uniforme sur [0,1], noté \mathcal{U}([0,1]) si sa fonction de répartition F_U est donnée par F_U(x) = \begin{cases} 0, & \text{si }x < 0\,, \\ x, & \text{si }x \in [0,1]\,, \\ 1, & \text{si }x > 1\,. \\ \end{cases}

Figure 1: Fonction de répartition de la loi uniforme

L’objectif est de simuler sur machine une suite U_1, \ldots, U_n de variables aléatoires i.i.d. de loi \mathcal{U}([0,1]). Plusieurs problèmes apparaissent alors :

  • Une machine est déterministe.
  • Les nombres entre 0 et 1 donnés par la machine sont de la forme k/2^p, pour k \in \{0, \ldots, 2^{p-1}\}. On ne pourra donc jamais générer des nombres qui ne sont pas de cette forme.
  • Vérifier qu’une suite est bien i.i.d. est un problème difficile.

Définition 1 (Générateur de nombres pseudo-aléatoires)

Un générateur de nombres pseudo-aléatoires (🇬🇧: Pseudo Random Number Generator, PRNG), est un algorithme déterministe récursif qui renvoie une suite U_1, \ldots, U_n dans [0,1] qui a un “comportement similaire” à une suite i.i.d. de loi \mathcal{U}([0,1]). Pour être plus rigoureux, ces nombres sont en fait des nombres entiers générés uniformément sur un certain interval. Dans un second temps, une transformation simple (normalisation) permet d’obtenir des nombres flottants (🇬🇧: floats) entre 0 et 1.

Pour aller plus loin

Parfois il est utile d’aller chercher dans le code source certaines information pour savoir comment les fonctions sont codées dans les packages que l’on utiliser. Par exemple, pour numpy que l’on utilise fréquement, on peut voir l’opération choisie ici: Random: int -> float en numpy.

Un tel algorithme se construit de la manière suivante :

  1. On part d’une graine (🇬🇧: seed) U_0 qui détermine la première valeur de manière la plus arbitraire possible.
  2. La procédure récursive s’écrit U_{n+1} = f(U_n), où f est une transformation déterministe, de sorte que U_{n+1} est le plus indépendant possible de U_1, \dots, U_n.

  • La fonction f est déterministe et prend ses valeurs dans un ensemble fini, donc l’algorithme est périodique. Le but est donc d’avoir la plus grande période possible.

  • Notons qu’une fois que la graine est fixée, alors l’algorithme donne toujours les mêmes valeurs. Fixer la graine peut donc être très utile pour répéter des simulations dans des conditions identiques et ainsi repérer des erreurs.

Exercice: bug ou feature?

Reprendre les widgets du chapitre Théorèmes asymptotiques et faites varier doucement le paramètre p (de Bernoulli). Que constatez-vous? Proposer une explication potentielle.

Générateur congruentiel linéaire

La plupart des PRNG s’appuient sur des résultats arithmétiques. Un des plus connus est celui appelé Générateur congruentiel linéaire (🇬🇧 Linear congruential generator, LCG). Il est défini comme suit: on construit récursivement une suite d’entiers X_i via la congruence X_{n+1} = a X_n + b \quad \text{mod } m \enspace, où a,b,m sont des entiers bien choisis pour que la suite obtenue ait de bonnes propriétés. Il suffit alors de considérer X_n/m. Par exemple, la fonction rand sur scilab utilise cette congruence avec m=2^{31}, a=843\; 314\; 861, et b=453\; 816\; 693.

Générateurs alternatifs

Les langages Python et R utilisent par défaut le générateur Mersenne-Twister qui s’appuie sur la multiplication vectorielle, mais d’autres générateurs sont aussi disponibles. Ce générateur a pour période m =2^{19937}-1, nombre qu’on peut raisonnablement considérer comme grand.

Pour numpy la méthode par défaut est PCG64 (cf. documentation de numpy), qui dispose de meilleures garanties statistiques (Voir le site https://www.pcg-random.org pour cela).

Usage en numpy

On suppose désormais disposer d’un générateur pseudo-aléatoire sur [0,1]. En numpy depuis la version 1.17, une bonne manière d’utiliser des éléments aléatoires est d’utiliser un générateur que l’on définit soi-même:

seed = 12345  # Toujours être conscient qu'une graine existe
rng = np.random.default_rng(seed)  #
print(rng.random())  ##  un tirage uniforme sur [0,1]
print(rng.random(size=5))  ## cinq tirages uniformes sur [0,1]
print(rng.random(size=(3, 2)))  ## matrice 3x2, à entrées unif. sur [0,1]
0.22733602246716966
[0.31675834 0.79736546 0.67625467 0.39110955 0.33281393]
[[0.59830875 0.18673419]
 [0.67275604 0.94180287]
 [0.24824571 0.94888115]]

Dans la suite on va voir comment générer d’autres lois à partir de la loi uniforme, mais il est clair que les logiciels modernes proposent un large éventail de distribution classique (gaussienne, exponentielle, etc.). Une liste exhaustive est donnée ici pour numpy.

Pour aller plus loin

Une excellent discussion sur les bonnes pratiques aléatoires en numpy, et l’usage de np.random.default_rng est donnée dans ce blog post d’Albert Thomas.

Propriété de la loi uniforme

On verra souvent apparaître la variable aléatoire 1-U où U \sim \mathcal{U}([0,1]). Il se trouve que 1-U suit aussi une loi uniforme sur [0,1] comme le montre le calcul de sa fonction de répartition. Ainsi pour tout x \in [0,1] on obtient

\begin{align*} \mathbb{P}(1-U \leq x) & = \mathbb{P}(U \geq 1-x),\\ & = 1-(1-x), \\ & = x\,. \end{align*} On peut démontrer facilement la même relation pour x<0 et x>1, d’où le résultat.

Méthode d’inversion

L’idée de la méthode d’inversion repose sur le résultat suivant :

Rappel sur la fonction quantile

Rappel : Pour F une fonction définie sur \mathbb{R} à valeurs dans [0, 1], croissante, on note

\forall q \in ]0,1[, \quad F^\leftarrow(q) = \inf\{ x \in \mathbb{R} : F(x)\geq q\} \tag{1}

Théorème 1 (Caractérisation des quantiles) Soit F une fonction définie sur \mathbb{R} à valeurs dans [0, 1], croissante et continue à droite, alors pour tout q \in ]0, 1[, on a

\begin{align} \{x \in \mathbb{R} : F(x) \geq q) \} & = \{x \in \mathbb{R} : x \geq F^\leftarrow(q) \} \end{align}

Preuve.

  • Cas \subset: Soit x \in \mathbb{R} t.q. F(x) \geq q, alors par définition de l’inf dans Équation 1, x \geq F^\leftarrow(q).
  • Cas \supset: Soit x \in \mathbb{R} t.q. x \geq F_X^\leftarrow(q) alors pour tout \epsilon > 0, x + \epsilon > F^\leftarrow(q), donc F(x + \epsilon) \geq q. Puis, par continuité à droite de F, F(x) \geq q.


Méthode d’inversion

Théorème 2 (Méthode d’inversion) Soit X une v.a réelle, et U \sim\mathcal{U}([0,1]), alors la variable aléatoire F_X^{\leftarrow}(U) a même loi que X.

Preuve. En utilisant le théorème précédent, on a \mathbb{P}(F_X^{\leftarrow}(U) \leq x) = \mathbb{P}(U \leq F_X(x)) pour tout x\in\mathbb{R}. Puis, comme U est une loi uniforme sur [0,1], \mathbb{P}(U\leq F_X(x))=F_X(x).

On en déduit donc que la loi de F_X^{-1}(U) est la même que celle de X, car les deux v.a. ont la même fonction de répartition.


Exemple 1 (Simulation d’une loi exponentielle) On rappelle que la loi exponentielle de paramètre \lambda > 0 a pour densité f_{\lambda}(x) = \lambda e^{-\lambda x} {1\hspace{-3.8pt} 1}_{\mathbb{R}_+}(x)\enspace. et donc pour fonction de répartition F_{\lambda}(x) = (1 - e^{-\lambda x}) {1\hspace{-3.8pt} 1}_{\mathbb{R}_+}(x)\enspace. On vérifie que F_{\lambda} est bijective de \mathbb{R}_+ dans ]0,1[ et que son inverse est donnée pour tout u \in ]0,1[ par F_{\lambda}^{-1}(u) = -\frac{1}{\lambda} \log(1-u)\enspace.


Exemple 2 (Simulation d’une loi de Weibull) La loi de Weibull de paramètre \lambda > 0 et k >0 est caractérisée par la fonction de répartition F(x) = (1 - e^{-(x/\lambda)^k}){1\hspace{-3.8pt} 1}_{\mathbb{R}_+}(x) C’est une loi utilisée dans différents domaines, notamment en gestion des risques (hydrologie, finance, assurance, etc.). Le calcul de l’inverse généralisée F^\leftarrow est immédiat car F est bijective sur ]0, \infty[ : F^\leftarrow (u) = \lambda (-\ln(1-u))^{\frac{1}{k}}\,, \quad u \in ]0,1[\,. La méthode d’inversion s’applique : si U \sim \mathcal{U}([0,1]), alors la variable aléatoire X = \lambda (-\ln(U))^{\frac{1}{k}} suit une loi de Weibull de paramètres \lambda et k.

Malheureusement, la fonction F n’est pas toujours inversible (penser aux lois discrètes) c’est donc pourquoi on utilise l’inverse l’inverse généralisée ou fonction quantile introduite dans la section Notations:

F^\leftarrow(p)= \inf\{ x \in \mathbb{R} \colon F(x)\geq p\} \enspace.

Interprétation: Définir l’inverse d’une fonction de répartition F revient à résoudre l’équation F(x) = \alpha d’inconnue x pour un \alpha fixé. Si F n’est pas bijective, deux problèmes apparaissent :

  • l’équation n’a aucune solution ce qui revient à dire que F n’est pas surjectif (graphiquement, F présente des sauts) ;
  • l’équation a plusieurs solutions ce qui revient à dire que F n’est pas injective (graphiquement cela se matérialise par un plateau à la hauteur \alpha). Un exemple classique est celui où F est la fonction de répartition d’une variable aléatoire discrète.

Le passage à l’inéquation F(x) \geq u permet de contourner la non-surjectivité : on ne regarde non plus les droites horizontales y=u mais la région \{y \geq \alpha\}. Le choix de l’\inf dans la définition de F^{\leftarrow} permet de contourner la non-injectivité : vu qu’il y a possiblement plusieurs x tels que F(x) \geq u, on choisit le “premier”. Ces considérations sont illustrées en Figure Figure 2.

Figure 2: Fonction de répartition et quantiles: solution de F(x) = u

Remarques additionnelles:

  • La fonction F étant croissante, la quantité F^\leftarrow(u) correspond au premier instant où F dépasse \alpha. Si F est bijective (ce qui équivaut dans ce cas à strictement croissante et injective), alors F^\leftarrow = F^{-1}.
  • La fonction F^\leftarrow n’est rien d’autre que la fonction quantile : si 0 < \alpha < 1, q_{1-\alpha} = F^\leftarrow(1-\alpha) est le quantile d’ordre (1-\alpha) de F. Par exemple, F^\leftarrow(1/2) correspond à la médiane.
  • Notons que si u=0, on peut alors naturellement poser F^{\leftarrow}(0) = -\infty. De même, avec la convention la convention \inf \emptyset = +\infty, on peut alors étendre la définition de F^\leftarrow à u=1 (mais F^\leftarrow(1) n’est pas toujours égal à \infty, voir les exemples ci-dessous).

Exemple 3 (Simulation d’une loi de Bernoulli) La fonction de répartition F d’une loi de Bernoulli de paramètre p \in ]0,1[ est donnée par F(x) = \Bigg\{ \begin{array}{ll} 0 & \text{ si } x < 0\,, \\ 1-p & \text{ si } 0 \leq x < 1\,, \\ 1 & \text{ si } x \leq 1\,. \end{array} L’inverse généralisée de F peut ainsi être calculée via la formule de l’exemple : F^\leftarrow(u) = \bigg\{ \begin{array}{ll} 0 & \text{ si } 0 < u \leq 1-p\,, \\ 1 & \text{ si } 1-p < u \leq 1\,. \end{array} ce qui se réécrit plus simplement F^\leftarrow(u) = {1\hspace{-3.8pt} 1}_{\{1-p < u\}}.

Ainsi, si U suit une loi uniforme sur [0,1] alors {1\hspace{-3.8pt} 1}_{\{1-p < U\}} suit une loi de Bernoulli de paramètre p. Comme 1-U suit aussi une loi uniforme sur [0,1], on en déduit que {1\hspace{-3.8pt} 1}_{\{U < p\}} suit une loi de Bernoulli de paramètre p. Notons que l’on peut remplacer l’inégalité stricte par une inégalité large.


Proposition 1 (Loi à support fini)

Soit X une variable aléatoire discrète prenant uniquement les valeurs x_1 < \dots < x_r (r modalité possibles) avec probabilité p_1, \dots, p_r (donc p_1 + \dots + p_r=1). On vérifie que pour tout u \in ]0,1[, F^\leftarrow(u) = \begin{cases} x_1 & \text{si } 0 < u \leq p_1\,, \\ x_2 & \text{si } p_1 < u \leq p_1+p_2\,, \\ & \vdots \\ x_r & \text{si } \sum_{i=1}^{r-1} p_i < u < 1\,. \end{cases}

Sur cet exemple, on peut prolonger la définition de F^\leftarrow à u=1 en posant F^\leftarrow(1) = x_r. L’inverse généralisée se réécrit alors sous la forme F^\leftarrow(u) = \sum_{k=1}^r x_k {1\hspace{-3.8pt} 1}_{ \{ \sum_{i=1}^{k-1}p_i < u \leq \sum_{i=1}^{k}p_i \} }\enspace, où on a posé p_0=0.

L’expression précédente s’étend directement au cas où X prend un nombre (infini) dénombrable de valeurs, la somme devenant alors une série.

La méthode est illustré ci-dessous pour quelques lois intéressantes:

Méthode de rejet

L’idée de la méthode de rejet est la suivante. On souhaite simuler une variable aléatoire X de densité f, appelée loi cible, mais f est trop compliquée pour que la simulation puisse se faire directement. On dispose cependant d’une autre densité g possédant les propriétés suivantes :

  • on sait simuler Y de loi g,
  • il existe m > 0 tel que f(x) \leq m \cdot g(x),
  • on sait évaluer le rapport d’acceptation r(x) = \frac{f(x)}{mg(x)}.

Remarquons d’ores et déjà que la constante m est nécessairement plus grande que 1 car 1 = \int_{\mathbb{R}} f(x) \, dx \leq m \int_{\mathbb{R}} g(x)\, dx = m\,.

L’idée est alors de considérer deux suites i.i.d. de variables aléatoires indépendantes entre elles:

  • (Y_n)_{n \geq 1} de loi g,
  • (U_n)_{n \geq 1} de loi uniforme sur [0,1].

En pratique, Y_n correspond à une proposition et U_n permettra de décider si on accepte la proposition ou non. Si oui, alors on conserve Y_n, sinon on simule Y_{n+1}. Le rapport d’acceptation, c’est-à-dire la proportion de Y_n acceptées, correspond à r(x).

Autrement dit, pour simuler X de densité f, il suffit de simuler Y de densité g et U uniforme jusqu’à ce que U \leq r(Y). La proposition suivante assure que cette méthode donne bien le résultat voulu.

Proposition 2 (Méthode de rejet)

Soit T = \inf \{n \geq 1 : U_n \leq r(Y_n)\} le premier instant où le tirage est accepté. Alors :

  • T suit une loi géométrique de paramètre 1/m,
  • la variable aléatoire X = Y_T a pour densité f et est indépendante de T.

Preuve. Il s’agit d’étudier la loi du couple (X,T). Pour x \in \mathbb{R} et n \in \mathbb{N}^{*}, on écrit \mathbb{P}(X \leq x, T=n)= \mathbb{P}(U_1 > r(Y_1), \dots, U_{n-1} > r(Y_{n-1}), U_n \leq r(Y_n), Y_n \leq x). Les n tirages étant iid, on obtient \mathbb{P}(X \leq x, T=n) = \mathbb{P}(U_1 > r(Y_1))^{n-1} \mathbb{P}(U_n \leq r(Y_n), Y_n \leq x)\,.

Concernant le premier terme, les variables aléatoires Y_1 et U_1 sont indépendantes donc leur loi jointe correspond au produit des densités : \begin{align*} \mathbb{P}(U_1 > r(Y_1)) & = \mathbb{P}((U_1, Y_1) \in \{(u,y) \in \mathbb{R}^2 : u > r(y)\}) \\ & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{u > r(y)\}} ({1\hspace{-3.8pt} 1}_{[0,1]}(u) g(y)) \, du dy \\ & = \int_{\mathbb{R}} \bigg( \int_0^1 {1\hspace{-3.8pt} 1}_{\{u > r(y)\}} \, du\bigg) g(y)\, d y \\ & = \int_{\mathbb{R}} (1-r(y)) \, g(y)\, d y\,, \end{align*} ce qui se réécrit, comme f et g sont des densités et que r(y) = f(y)/(m \cdot g(y)): \begin{align*} \mathbb{P}(U_1 > r(Y_1)) & = \int_{\mathbb{R}} g(y)\, d y - \int_{\mathbb{R}} \dfrac{f(y)}{m}\, dy \\ & = 1 - \dfrac{1}{m}\,. \end{align*} Le deuxième terme se calcule de manière analogue :

\begin{align*} \mathbb{P}(U_n \leq r(Y_n), Y_n \leq x) & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{u \leq r(y)\}} {1\hspace{-3.8pt} 1}_{\{y \leq x\}} ({1\hspace{-3.8pt} 1}_{[0,1]}(u) g(y)) \, du dy \\ & = \int_{\mathbb{R}} \bigg( \int_0^1 {1\hspace{-3.8pt} 1}_{\{u \leq r(y)\}} \, du\bigg) {1\hspace{-3.8pt} 1}_{\{y \leq x\}} g(y)\, d y\,, \end{align*} c’est-à-dire

\begin{align*} \mathbb{P}(U_n \leq r(Y_n), Y_n \leq x) & = \int_{\mathbb{R}} r(y) {1\hspace{-3.8pt} 1}_{\{y \leq x\}} g(y)\, d y \\ & = \int_{-\infty}^x \dfrac{f(y)}{m}\, d y \\ & = \dfrac{F(x)}{m}\,, \end{align*} où F est la fonction de répartition de la loi de densité f. On peut ainsi conclure que \mathbb{P}(X \leq x, T=n) = \bigg(1 - \dfrac{1}{m}\bigg)^{n-1} \dfrac{F(x)}{m}\,. Il ne reste plus qu’à étudier les lois marginales. D’une part, par continuité monotone croissante, \mathbb{P}(T=n) = \lim_{q \to \infty} \mathbb{P}(X \in ]-\infty, q], T=n)\,, ce qui donne

\begin{align*} \mathbb{P}(T=n) & = \lim_{q \to \infty} \bigg(1 - \dfrac{1}{m}\bigg)^{n-1} \dfrac{F(q)}{m}\\ & = \bigg(1 - \dfrac{1}{m}\bigg)^{n-1} \dfrac{1}{m}\,. \end{align*} On en déduit que T suit une loi géométrique de paramètre 1/m. D’autre part, par \sigma-additivité,

\begin{align*} \mathbb{P}(X \leq x) & = \mathbb{P}(X \leq x, T \in \mathbb{N}^*)\\ & = \sum_{n=1}^\infty \mathbb{P}(X \leq x, T=n)\,, \end{align*} ce qui donne

\begin{align*} \mathbb{P}(X \leq x) & = \sum_{n=1}^\infty \bigg(1 - \dfrac{1}{m}\bigg)^{n-1} \dfrac{F(x)}{m}\\ & = \dfrac{1}{1-(1-1/m)} \dfrac{F(x)}{m}\\ & = F(x)\,, \end{align*} ce qui prouve que X a pour loi F.

Enfin, la loi du couple (X,T) est égale au produit des lois

\begin{align*} \mathbb{P}(X \leq x, T=n) & = \bigg(1 - \dfrac{1}{m}\bigg)^{n-1} \dfrac{F(x)}{m}\\ & = \mathbb{P}(T=n) \mathbb{P}(X \leq x)\,, \end{align*}

def accept_reject(n, f, g, g_sampler, m):
    """
    n: nombre de simulations
    f: densité cible
    g: densité des propositions, g_sampler: simulateur selon g
    m: constante pour la majoration
    """
    x_samples = np.zeros(n)
    u_samples = np.zeros(n)
    accepted = np.zeros(n)
    n_accepted = 0
    while n_accepted < n:
        x = g_sampler()
        u = np.random.uniform()
        alpha = u * m * g(x)
        u_samples [n_accepted] = alpha
        x_samples[n_accepted] = x
        if  alpha <= f(x):
            accepted[n_accepted] = 1
        n_accepted += 1
    return x_samples, u_samples, accepted
En pratique…

On simule U_1 et Y_1. Si U_1 \leq r(Y_1) c’est gagné, on pose X=Y_1. Sinon, on simule U_2 et Y_2 et on teste à nouveau l’inégalité U_2 \leq r(Y_2). Et ainsi de suite. Comme T suit une loi géométrique de paramètre 1/m, son espérance vaut m : il faut en moyenne m tentatives pour obtenir une simulation de la loi de densité f. L’objectif est alors de choisir un couple (g, m) de sorte que m soit le plus proche possible de 1.

Exemple 4 (Rejet d’une loi polynomiale) Donnons un exemple jouet (on étudiera des exemples plus pertinents en TD). On considère la densité f(x) = 4x^3 {1\hspace{-3.8pt} 1}_{[0,1]}(x). Comme f est majorée par 4, on peut choisir pour g la densité de la loi uniforme sur [0,1] et m=4. Alors, r(x) =f(x) / (mg(x)) = x^3, pour x \in [0,1]. On simule donc (Y_1, U_1) et on teste si U_1 \leq Y_1^3, etc.

Bien évidemment, on privilégiera ici une simulation via F^\leftarrow qui permet de générer des variables aléatoires de loi f plus rapidement.

Figure 3: Visualisation des zones d’acceptations/rejet (g uniforme)

Nous pouvons facilement améliorer la proportion de point acceptés en proposant par exemple g définie par g(x) = 2x {1\hspace{-3.8pt} 1}_{[0, 1]}(x), et m=2.

Figure 4: Visualisation des zones d’acceptations/rejet (g triangulaire)


Exemple 5 (Rejet d’une loi de densité d’Andrews) Considérons la densité d’Andrews définie par f(x) = \frac{1}{S} \frac{\sin(\pi\cdot x)}{\pi \cdot x} {1\hspace{-3.8pt} 1}_{[-1,1]}(x), avec S = \int_{-1}^{1}\frac{\sin(\pi\cdot x)}{\pi \cdot x}dx. Dans ce contexte, on ne connait pas la valeur exacte de S, et on va donc utiliser la méthode de rejet pour simuler des variables aléatoires de loi f sans cette information. On peut l’adapter le test de la manière suivante: si l’on prend m=2/S et g(x) = \frac{1}{2} {1\hspace{-3.8pt} 1}_{[-1,1]}(x), on observe que tester u\leq \frac{f(x)}{m \cdot g(x)} est équivalent à tester u \leq r(x)=\frac{1}{S} \frac{\sin(\pi\cdot x)}{\pi \cdot x} \cdot \frac{1}{\frac{2}{S} g(x)} = \frac{\sin(\pi\cdot x)}{\pi \cdot x} \cdot \frac{1}{2 \cdot g(x)}, ce qui peut se faire sans connaissance de S. De plus on peut vérifier que g(x) = \frac{1}{2} {1\hspace{-3.8pt} 1}_{[-1,1]}(x) définit une densité et que f(x) \leq m \cdot g(x) pour tout x\in \mathbb{R}.

n = 10000
g = lambda x: np.ones_like(x) / 2
g_sampler = lambda: 2 * np.random.uniform() - 1
m = 2

x_samples, u_samples, accepted = accept_reject(n, np.sinc, g, g_sampler, m)
ratio = np.sum(accepted) / n
# Note: https://stackoverflow.com/questions/70804891/how-to-vectorize-a-function-in-python-that-includes-a-limit-e-g-sinx-x

On peut approcher numériquement la valeur exacte de S en utilisant une méthode de calcul approchée, ce qui permet de comparer ici notre méthode de rejet avec la densité sous-jacente:

from scipy import integrate
S = integrate.quad(np.sinc, -1, 1)[0]
print(f"En utilisant la méthode de rejet, on trouve que S = {S:.3f}")
En utilisant la méthode de rejet, on trouve que S = 1.179

Enfin, on peut visualiser la qualité l’approximation de la densité par la méthode de rejet en comparant la densité approchée (avec un histogramme) avec la densité exacte:

Figure 5: Méthode de rejet pour simuler une loi de densité de type Andrews, sans connaissance de la valeur exacte de la constante de normalisation.

Cas mutlidimensionnel

Commençons par un cas de dimension deux.

Pour cela on va utiliser la méthode de rejet pour simuler une loi de densité f sur \mathbb{R}^2. En particulier, un exemple classique est de tirer des points dans le disque unité, c’est-à-dire de simuler une loi uniforme sur le disque unité. Pour cela, on va utiliser la méthode de rejet avec g la densité de la loi uniforme sur le carré [-1,1]^2.

Mais prenons un autre exemple, à savoir tirer des points uniformément dans la surface délimité par une cardioïde. Pour cela, on va utiliser la méthode de rejet avec g la densité de la loi uniforme sur le carré [-2,2]^2.

Figure 6: Méthode de rejet pour simuler une loi uniforme sur un disque unité. 
Aire estimée: 3.1493333333333333
Figure 7: Méthode de rejet pour simuler une loi uniforme sur une surface délimitée par une cardioïde. 
Aire estimée: 4.890000000000001
EXERCICE loi uniforme sur un cylindre

Proposer une méthode pour simuler une loi uniforme sur un cylindre de rayon 1 et de hauteur 10.

Autres méthodes

Sommation de variables aléatoires

Pour simuler une variable aléatoire de loi binomiale \mathcal{B}(n,p), on peut utiliser la méthode d’inversion. Cependant, cela nécessite le calcul de l’inverse généralisée de F, donc de coefficients binomiaux et de puissances de p et 1-p. À la place, on utilisera plutôt la relation bien connue suivante : si X_1, \ldots, X_n est une suite iid de variables aléatoires de loi de Bernoulli de paramètre p, alors X = X_1 + \cdots + X_n \sim \mathcal{B}(n,p)\,.

Pour simuler des variables aléatoires de Bernoulli, on utilise la méthode d’inversion (voir Exemple ). Ainsi, si U_1, \ldots, U_n sont des variables aléatoires iid de loi uniforme sur [0,1], alors \sum_{i=1}^n {1\hspace{-3.8pt} 1}_{\{U_i \leq p\}} \sim \mathcal{B}(n,p)\,.

Loi de Poisson

Rappelons qu’une variable aléatoire X suit une loi de Poisson de paramètre \lambda > 0, notée X \sim \mathcal{P}(\lambda) si \mathbb{P}(X = k) = e^{-\lambda} \dfrac{\lambda^k}{k!}\,, \quad k \in \mathbb{N}\,. Une méthode pour simuler une variable aléatoire de loi de Poisson est donnée par la proposition suivante.

Proposition 3 (Génération de v.a. de loi de Poisson)

Soit (E_n)_{n \geq 1} des variables aléatoires i.i.d. de loi exponentielle de paramètre \lambda > 0. On pose S_k = E_1 + \cdots + E_k. Alors, pour tout n \in \mathbb{N} \mathbb{P}(S_n \leq 1 < S_{n+1}) = e^{-\lambda} \dfrac{\lambda^n}{n!}\enspace . Ainsi, la variable aléatoire T définie par T \triangleq \sup \{n \in \mathbb{N} : S_n \leq 1\} suit une loi de Poisson de paramètre \lambda : T \sim \mathcal{P}(\lambda).

La preuve repose sur le lemme suivant.

Lemme 1 (Loi de Erlang)

Soit n variables aléatoires E_1, \dots, E_n i.i.d. de loi exponentielle de paramètre \lambda >0. La somme E_1+\dots+E_n suit une loi d’Erlang de paramètres (n,\lambda), donnée par la fonction de répartition F_{n,\lambda}(t) = 1 - \sum_{k=0}^{n-1} e^{-\lambda t} \frac{(\lambda t)^k}{k!}\,.

Preuve. On montre le résultat pour n=2. La généralisation à k quelconque se fait par récurrence. Soit t > 0, et f_{\lambda}(x)={1\hspace{-3.8pt} 1}_{\{x \geq 0 \}} \lambda e^{-\lambda x} la densité d’une loi exponentielle de paramètre \lambda. Les variables aléatoires E_1 et E_2 étant indépendantes et suivant des lois exponentielles de paramètre \lambda_1 et \lambda_2, on a

\begin{align*} \mathbb{P}(E_1+E_2 \leq t) & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{x_1 + x_2 \leq t\}} f_{\lambda}(x_1) f_{\lambda}(x_2)\, d x_1 d x_2 \\ & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{x_1 + x_2 \leq t\}} \lambda^2 e^{-\lambda (x_1+x_2)} {1\hspace{-3.8pt} 1}_{\{x_1 \geq 0\}} {1\hspace{-3.8pt} 1}_{\{x_2 \geq 0\}}\, d x_1 d x_2 \\ & = \int_{\mathbb{R}^2} {1\hspace{-3.8pt} 1}_{\{0 \leq x_1 \leq t\}} {1\hspace{-3.8pt} 1}_{\{0 \leq x_2 \leq t-x_1\}} \lambda^2 e^{-\lambda x_1} e^{-\lambda x_2}\, d x_1 d x_2 \\ & = \int_0^t \lambda e^{-\lambda x_1} \bigg(\int _0^{t-x_1} \lambda e^{-\lambda x_2}\, d x_2\bigg) d x_1\,. \end{align*} La première intégrale se calcule alors facilement : \int _0^{t-x_1} \lambda e^{-\lambda x_2}\, d x_2 = 1 - e^{-\lambda(t-x_1)}\,. On obtient alors

\begin{align*} \mathbb{P}(E_1+E_2 \leq t) & = \int_0^t \lambda e^{-\lambda x_1}dx_1 - \int_0^t e^{-\lambda t} d x_1\\ & = 1 - e^{-\lambda t} - \lambda t e^{-\lambda t}\,. \end{align*} Si t<0, alors comme les E_i ne prennent que des valeurs positives on trouve \mathbb{P}(E_1 + E_2 \leq t) = 0. Ceci prouve le résultat pour n=2.

On peut désormais prouver le résultat de la Proposition 3.

Preuve. Pour n \in \mathbb{N}, on décompose la probabilité \mathbb{P}(S_n \leq 1 < S_{n+1}) via

\begin{align*} \mathbb{P}(S_n \leq 1 < S_{n+1}) & = \mathbb{P}(\{S_n \leq 1\} \setminus \{S_{n+1} \leq 1\})\\ & = \mathbb{P}(S_n \leq 1) - \mathbb{P}(S_{n+1} \leq 1)\,. \end{align*} Le lemme précédent donne \mathbb{P}(S_n \leq 1) = 1 - \sum_{k=0}^{n-1} e^{-\lambda} \dfrac{\lambda^k}{k!} et \mathbb{P}(S_{n+1} \leq 1) = 1 - \sum_{k=0}^{n} e^{-\lambda} \dfrac{\lambda^k}{k!}\,. On obtient alors le résultat souhaité : \mathbb{P}(S_n \leq 1 < S_{n+1}) = e^{-\lambda} \dfrac{\lambda^n}{n!}\,.

On conclut la preuve de la proposition en remarquant que \mathbb{P}(T=n) = \mathbb{P}(S_n \leq 1 < S_{n+1})\,.

La simulation d’une variable aléatoire de Poisson repose donc sur la simulation de lois exponentielles qui se fait via la méthode d’inversion, comme vu dans Exemple 1. En pratique, on simule E_1 et on teste si E_1 > 1. Si oui, on pose alors T=0. Si non, on simule E_2 et on teste si E_1 + E_2 > 1. Si oui, on pose T=1. Sinon on continue la procédure.

Bibliographie et pour aller plus loin

Retour au sommet