15 minutes de Python

Méthode d’inversion

La méthode d’inversion permet de simuler une variable aléatoire réelle à densité à partir d’une variable uniforme. Dans toute cette page, on se place dans le cadre suivant :

• \(X\) est une variable aléatoire réelle admettant une densité \(f\) et de fonction de répartition \( F\),
• la densité \( f \) est nulle en dehors d’un intervalle \(\left]a,b\right[\) (avec \(a\) et \(b\) éventuellement infinis),
• la fonction de répartition \(F\) est strictement croissante sur \(\left]a,b\right[\).

Fondements mathématiques

Bijection induite par la fonction de répartition

Sur l’intervalle \(\left]a,b\right[\), la fonction de répartition \(F\) est continue et strictement croissante. De plus, \[ \lim_{x\to a^+}F(x)=0 \qquad\text{et}\qquad \lim_{x\to b^-}F(x)=1 \]

Ainsi, l’application \[ G : \left]a,b \right[ \to \left]0,1 \right[, \quad x \mapsto F(x) \] est une bijection. On peut donc définir sa fonction réciproque \[ G^{-1} : \left]0,1\right[ \to \left]a,b\right[ \]

Loi de \(G^{-1}(U)\)

Soit \(U\sim\mathcal{U}(\left]0,1\right[)\) et posons \[ Y = G^{-1}(U) \] Pour tout \(x\in\left]a,b\right[\), comme \(G\) est strictement croissante, on a \begin{align*} Y \le x &\Leftrightarrow\; G^{-1}(U) \le x \\ &\Leftrightarrow\; U \le G(x) \end{align*}

Donc, comme \( G(x) = F(x) \) lorsque \( x \in \left] a,b \right[ \) : \begin{align*} \mathbb{P}(Y\le x) &=\mathbb{P}(U\le G(x)) \\ &=\mathbb{P}(U\le F(x)) \\ &=F(x) \end{align*} car \(U\) est uniforme sur \(\left]0,1\right[\) et \( F(x) \) appartient à \( ]0,1[ \). Ainsi, \(Y\) admet la même fonction de répartition que \(X\), donc elle suit la même loi que \( X\).

Conséquence pour la simulation

Le raisonnement précédent montre que si l’on sait calculer explicitement la fonction réciproque \(G^{-1}\), alors on peut simuler \(X\) de la façon suivante :

1. on simule une variable uniforme \(U\) sur \(\left]0,1\right[\),
2. on calcule \(X = G^{-1}(U)\).

C’est exactement le principe de la méthode d’inversion.

Exemple : loi exponentielle

Soit \(X\) une variable aléatoire de loi exponentielle de paramètre \(\lambda>0\). Elle admet une densité nulle en dehors de \(]0,+\infty[\), et sa fonction de répartition vérifie \[ F(x)=\left\{ \begin{array}{ll} 0 & \text{si }x\le 0\\ 1-e^{-\lambda x} & \text{si }x>0 \end{array} \right. \] Sur \(]0,+\infty[\), \(F\) est strictement croissante.

Pour \(u\in\left]0,1\right[\) et \( x \in \left] 0, + \infty \right[ \), on a de plus : \begin{align*} u=1-e^{-\lambda x} &\Leftrightarrow e^{-\lambda x} = 1-u &\\ &\Leftrightarrow x=-\dfrac{1}{\lambda} \, \ln(1-u) \end{align*} Ainsi, si \(U\sim\mathcal{U}(\left]0,1\right[)\), alors \[ X=-\dfrac{1}{\lambda} \, \ln(1-U) \] suit une loi exponentielle de paramètre \(\lambda\) d’après le point précédent et on peut donc simuler une réalisation de \( X\) avec le code suivant :

import numpy as np
import numpy.random as rd

# Paramètre
Lambda = 2.0

# 1) Simulation d'une uniforme sur ]0,1[
U = rd.random(N)

# 2) Inversion : X = G^{-1}(U)
X = -np.log(1 - U) / Lambda

# Renvoie de la valeur
print(X)

import numpy as np
import numpy.random as rd

alpha = 3

U = rd.random()

X = 1 / (1 - U)**(1/alpha)

print("U =", U)
print("X =", X)