seed = 12345 # choix de la graine
rng = np.random.default_rng(seed) # générateurHAX603X: Modélisation stochastique
2024-01-19
Question: Comment simuler en pratique des variables aléatoires i.i.d?
Approche: Commencer par les v.a. uniformes et en déduire les autres lois
Rappel : \(U\) suit une loi uniforme sur \([0,1]\): \(U\sim\mathcal{U}([0,1])\) ssi sa fonction de répartition \(F_U\) vaut \[ F_U(x) = \begin{cases} 0, & \text{si }x < 0\,, \\ x, & \text{si }x \in [0,1]\,, \\ 1, & \text{si }x > 1\,. \end{cases} \]
Objectif: simuler sur machine une suite \(U_1, \dots, U_n\) de v.a., i.i.d., de loi \(\mathcal{U}([0,1])\).
Difficultés:
Définition 1 (Générateur de nombres pseudo-aléatoires)
Un générateur de nombres pseudo-aléatoires (:gb:: 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])\).
Remarque: ces nombres sont obtenus depuis des nombres entiers générés aléatoirement et uniformément sur grand interval, puis une transformation simple (normalisation) permet d’obtenir des nombres flottants (:gb:: floats) entre 0 et 1.
Un PRNG se construit ainsi :
La plupart des PRNG s’appuient sur des résultats arithmétiques.
Le plus célèbre: Générateur Congruentiel Linéaire (:gb: Linear Congruential Generator, LCG).
Récurrence: \[ X_{n+1} = a X_n + b \quad \text{mod } m \enspace, \] \(a,b,m\), entiers bien choisis pour que la suite obtenue ait de bonnes propriétés
Normalisation: \(X_n/m\).
Exemple: la fonction rand de scilab utilisait cette congruence avec \(m=2^{31}\), \(a=843\; 314\; 861\), et \(b=453\; 816\; 693\).
Python et R: Mersenne-Twister, s’appuie sur la multiplication vectorielle (période du générateur \(m =2^{19937}-1\))numpy: PCG64 (cf. documentation de numpy), dispose de meilleures garanties statistiques; voir https://www.pcg-random.orgnumpyEn numpy (version>1.17): utiliser des éléments aléatoires est d’utiliser un générateur
[0.31675834 0.79736546 0.67625467 0.39110955 0.33281393]Suite du cours: apprendre à générer de nombreuses lois à partir de la loi uniforme
En pratique: les logiciels proposent les distributions classiques (gaussiennes, exponentielles, etc.), utiliser plutôt ces fonctions que de les implémenter soi-même.
Liste exhaustive pour numpy:
https://numpy.org/doc/stable/reference/random/generator.html#distributions
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.
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}\]
Note: Si \(x_0 \in \{ x \in \mathbb{R} : F(x)\geq q\}\) alors \([x_0,+\infty[ \subset \{ x \in \mathbb{R} : F(x)\geq q\}\)
Théorème 1 (Caraté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} \]
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\): Soient \(\epsilon>0\) et \(x \in \mathbb{R}\) t.q. \(x \geq F^\leftarrow(q)\) alors (def. de l’inf) \(\exists x_0 \in \{ x \in \mathbb{R} : F(x)\geq q\}\), t.q. \(x + \epsilon > x_0\). Ainsi, \(F(x + \epsilon) \geq F(x_0) \geq q\); par continuité à droite de \(F\), \(F(x) \geq q\)
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 \[ \forall x\in\mathbb{R}, \quad \mathbb{P}(x \geq F_X^\leftarrow(U)) = \mathbb{P}(F_X(x) \geq U) \]
Puis, comme \(U\) est une loi uniforme sur \([0,1]\):
\[ \mathbb{P}(F_X(x) \geq U) = F_X(x) \]
Ainsi, \(F_X^\leftarrow(U)\) et \(X\) ont même loi: les deux v.a. ont la même fonction de répartition
Proposition 1 (Symétrie de la loi uniforme) Soit \(U \sim \mathcal{U}([0,1])\) une variable uniforme sur \([0,1]\). Alors, \(1-U\) suit aussi une loi uniforme sur \([0,1]\).
Preuve:
On va décrire la fonction de répartition de \(1-U\) et montrer qu’elle est égale à celle d’une loi uniforme sur \([0,1]\).
Le résultat est facile pour \(x \notin [0,1]\), on suppose donc \(x \in [0,1]\).
\[ \begin{align*} \mathbb{P}(1-U \leq x) & \class{fragment}{{}= \mathbb{P}(U \geq 1-x) }\\ & \class{fragment}{{} = 1-(1-x)} \\ & \class{fragment}{{} = x} \end{align*} \]
Avec le résultat: \[ U \sim \mathcal{U}([0,1]) \iff 1-U \sim \mathcal{U}([0,1])\enspace, \]
Pour simuler une loi exponentielle: simuler \(U\) uniforme et appliquer \(-\tfrac{1}{\lambda} \log(\cdot)\)
\[ \boxed{-\tfrac{1}{\lambda} \log(U) \sim \mathcal{E}(\lambda)} \]
Voir animation dans la section Cours, section “Méthode d’inversion”.
Motivation: simuler une variable aléatoire \(X\) de densité \(f\) (loi cible), mais \(f\) est trop compliquée pour la méthode de l’inverse.
Idée: tirer suivant une autre loi \(g\) (loi des propositions) et rejeter certains tirages.
Remarque 1: \(g(x)=0 \implies f(x)=0\), ainsi le support de \(g\) doit englober celui de \(f\)
Remarque 2: \(m \geq 1\) car \(m = m \displaystyle\int_\mathbb{R} g(x)\, dx \class{fragment}{{} \geq \displaystyle\int_\mathbb{R} f(x) dx} \class{fragment}{{} = 1}\)
Considérer deux suites i.i.d. de v.a. indépendantes entre elles:
En pratique, \(Y_n\) correspond à une proposition et \(U_n\) permettra de décider l’acceptation/rejet de la proposition:
Pour simuler \(X\) de densité \(f\), simuler \(Y_n\) (suivant \(g\)), \(U_n\) (suivant \(\mathcal{U}[0,1]\)) et accepter si \[ U_n \leq r(Y_n) = \frac{f(Y_n)}{m\cdot g(Y_n)} \]
Note
Dans la suite on verra qu’on tire des points \((Y_n, 4U_n)\) et qu’on teste si ils sont dans l’ensemble \(\{(x,y) \in \mathbb{R}^2: y \leq f(x) \}\)
def accept_reject(n, f, g, g_sampler, m, rng):
"""
n: nombre de simulations
f: loi cible
g, g_sampler: loi et générateur des propositions
m: constante pour la majoration
rng: générateur pseudo-aléatoire
"""
x_samples = np.zeros(n)
u_samples = np.zeros(n)
accepted = np.zeros(n)
for i in range(n):
x = g_sampler()
u = rng.uniform()
alpha = u * m * g(x) # note: pour le test on peut éviter les divisions
u_samples [i] = alpha
x_samples[i] = x
if alpha <= f(x):
accepted[i] = 1
return x_samples, u_samples, accepted
Supposons disposer d’un générateur de loi triangulaire sur \([0,1]\) (cf. np.random.triangular(0, 1, 1))
Conclusion: plus le majorant est proche de la loi cible, plus le taux d’acceptation est élevé, et moins de simulations sont nécessaires
Note
L’exemple est pour l’illustration de la méthode, dans le cas présent méthode de l’inverse fonctionnerait aussi (on peut calculer la fonction quantile explicitement).
Rappel:
Théorème 3 (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 :
Pour \(x \in \mathbb{R}\) et \(n \in \mathbb{N}^{*}\), et \(X = Y_T\), 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)\)
\[ \mathbb{P}(X \leq x, T=n) = {\color{blue}\mathbb{P}(U_1 > r(Y_1))}^{n-1} \cdot {\color{brown}\mathbb{P}(U_n \leq r(Y_n), Y_n \leq x)}, \quad \textsf{( tirages i.i.d.)} \]
Premier terme: \(Y_1\) et \(U_1\) sont indépendantes, leur loi jointe correspond au produit des densités :
\[ {\color{blue} \begin{align*} \mathbb{P}(U_1 > r(Y_1)) & \class{fragment}{{} = \mathbb{P}((U_1, Y_1) \in \{(u,y) \in \mathbb{R}^2 : u > r(y)\})} \\ & \class{fragment}{{} = \int_{\mathbb{R}^2} \left( {1\hspace{-3.8pt} 1}_{\{u > r(y)\}} \right) \cdot \left({1\hspace{-3.8pt} 1}_{[0,1]}(u) g(y)\right) \, du dy} \\ & \class{fragment}{{} = \int_\mathbb{R} \left( \int_0^1 {1\hspace{-3.8pt} 1}_{\{u > r(y)\}} \, du\right) g(y)\, d y} \class{fragment}{{} = \int_\mathbb{R} (1-r(y)) \, g(y)\, d y}\\ & \class{fragment}{{} = \int_\mathbb{R} g(y) - \int_\mathbb{R}\frac{f(y)}{m} d y, \quad\quad \text{car } r(y) = \frac{f(y)}{m \cdot g(y)}}\\ & \class{fragment}{{} = 1-\tfrac{1}{m}} \end{align*} } \]
Second terme: \[ {\color{brown} \begin{align*} \mathbb{P}(U_n \leq r(Y_n), Y_n \leq x) & \class{fragment}{{ = \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 }} \\ & \class{fragment}{{} = \int_\mathbb{R} \left( \int_0^1 {1\hspace{-3.8pt} 1}_{\{u \leq r(y)\}} \, du\right) {1\hspace{-3.8pt} 1}_{\{y \leq x\}} g(y)\, d y} \\ & \class{fragment}{{} = \int_\mathbb{R} r(y) {1\hspace{-3.8pt} 1}_{\{y \leq x\}} g(y)\, d y } \\ & \class{fragment}{{} = \int_{-\infty}^x \dfrac{f(y)}{m}\, d y } \\ & \class{fragment}{{} = \dfrac{F(x)}{m} , \quad F \textsf{ fonction de répartition associée à} f} \end{align*} } \]
\[ \mathbb{P}(X \leq x, T=n) = {\color{blue}\left(1 - \tfrac{1}{m}\right)^{n-1}} \cdot {\color{brown}\tfrac{F(x)}{m}} \]
On peut alors obtenir les lois marginales: \[ \begin{align*} \mathbb{P}(T=n) & = \lim_{q \to \infty} \mathbb{P}(X \in ]-\infty, q], T=n) = \lim_{q \to \infty} \left(1 - \tfrac{1}{m}\right)^{n-1} \tfrac{F(q)}{m}\\ & = \left(1 - \tfrac{1}{m}\right)^{n-1} \tfrac{1}{m} \end{align*} \] Ainsi, \(T\) suit une loi géométrique de paramètre \(1/m\), puis \(X\) a pour loi \(F\):
\[ \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) \\ & = \sum_{n=1}^\infty \left(1 - \tfrac{1}{m}\right)^{n-1} \tfrac{F(x)}{m} = \tfrac{1}{1-(1-1/m)} \tfrac{F(x)}{m} = F(x) \enspace. \end{align*} \]
On obtient l’indépendance de \(T\) et \(X\) car on peut alors écrire: \[ \forall x \in \mathbb{R}, \forall n \in \mathbb{N}^*, \quad \mathbb{P}(X \leq x, T=n) = \mathbb{P}(X \leq x) \cdot \mathbb{P}(T=n) \]
Loi de Andrews (densité proportionnelle à \(\mathrm{sinc}\), sinus cardinal): \[ \forall x \in \mathbb{R},\quad 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\) non explicite. On note parfois: \(f(x) \propto \frac{\sin(\pi\cdot x)}{\pi \cdot x} {1\hspace{-3.8pt} 1}_{[-1,1]}(x)\)
Méthode du rejet: prendre \(m=2/S\) et \(g(x) = \frac{1}{2} {1\hspace{-3.8pt} 1}_{[-1,1]}(x)\):
\[ u \leq r(x) = \frac{f(x)}{m \cdot g(x)} \iff u \leq \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)} \]
Ainsi l’évaluation de \(r(x)\) est possible sans connaître \(S\)!
Soit \(\tilde{f}: \mathbb{R} \to [0,+\infty[\) connue et \(S \triangleq \int_{\mathbb{R}} \tilde{f}(x) \, d x < + \infty\) inconnue (ou dure à évaluer)
Densité cible: \(\quad f(x) = \frac{\tilde{f}(x)}{S}\)
Méthode du rejet pour \(f\), en utilisant seulement \(\tilde{f}\): soit \(\tilde{m}>0\) un majorant de \(\tilde{f}\) t.q. \[ \begin{align*} \tilde{f}(x) \leq \tilde{m} \cdot g(x) \end{align*} \]
Application avec \(m=\tilde{m}/S\) (sans connaître \(S\)), le test d’acceptation donne:
\[ \begin{align*} U_n & \leq \frac{f(Y_n)}{m \cdot g(Y_n)}\\ & \class{fragment}{{}\leq \frac{\frac{\tilde{f}(Y_n)}{S}}{\frac{\tilde{m}}{S} \cdot g(Y_n)}} \class{fragment}{{}= \frac{\tilde{f}(Y_n)}{\tilde{m} \cdot g(Y_n)}} \end{align*} \]
impossibilité de la méthode de l’inverse: fonction de répartition non disponible (en général)
la méthode de rejet : généralisable au cas multidimensionnel
Loi cible: loi uniforme sur le disque unité, \(f(x)\propto {1\hspace{-3.8pt} 1}_{x_1^2+x_2^2 \leq 1}(x)\) pour \(x=(x_1,x_2)\in\mathbb{R}^2\)
Loi majorante: loi uniforme sur le carré \([-1,1]^2\), \(g(x)\triangleq \tfrac{1}{2} {1\hspace{-3.8pt} 1}_{[-1, 1]^2}(x)\) et \(m=2\)
Note
La loi uniforme sur le carré est une loi produit: il suffit de savoir générer une loi uniforme sur un segment 1D pour l’obtenir
Loi cible: loi uniforme sur le disque unité, \(f(x)\triangleq{1\hspace{-3.8pt} 1}_{(x_1^2+x_2^2 - x_2)^2 \leq x_1^2+ x_2^2}(x)\) pour \(x=(x_1,x_2)\in\mathbb{R}^2\)
Loi majorante: loi uniforme sur le rectangle \([-2,3]\times [-1.5,1.5]\), \(g(x)\triangleq \tfrac{1}{15}{1\hspace{-3.8pt} 1}_{[-2,3]\times [-1.5,1.5]}(x)\) et \(m=15\)
Note
La loi uniforme sur un rectangle est une loi produit: il suffit de savoir générer une loi uniforme sur un segment 1D pour l’obtenir
Théorème 4 (Géneration uniforme sur un ensemble) Soient \(A\subset B \subset \mathbb{R}^d\), deux ensembles mesurables pour la mesure de Lebesgue. Pour générer selon \(\mathcal{U}(A)\), connaissant un générateur selon \(\mathcal{U}(B)\), la méthode du rejet consiste ici à tirer \(Y_i \sim \mathcal{U}(B)\) (i.i.d) et à garder \(Y_i\) si \(Y_i \in A\).
Preuve: On note \(f\triangleq \frac{1}{|A|} {1\hspace{-3.8pt} 1}_{A}\), \(g\triangleq \frac{1}{|B|}{1\hspace{-3.8pt} 1}_{B}\) et \(m\triangleq\frac{|B|}{|A|}\). Comme \(A \subset B\), pour tout \(x\in \mathbb{R}^{d}\): \[ \begin{align*} \class{fragment}{{}f(x)} \class{fragment}{{}= \frac{1}{|A|} {1\hspace{-3.8pt} 1}_{A}(x)} \class{fragment}{{}\leq \frac{1}{|B|} {1\hspace{-3.8pt} 1}_{B}(x) \cdot \frac{|B|}{|A|}} \class{fragment}{{}\leq g(x) \cdot m} \end{align*} \]
\[ \begin{align*} &\class{fragment}{{}Y \sim \mathcal{U}(B),} \quad \class{fragment}{{} U \sim \mathcal{U}([0,1])} \\ &\class{fragment}{{}r(Y)} \class{fragment}{{}= \frac{f(Y)}{m\cdot g(Y)}} \class{fragment}{{}= \frac{ \frac{1}{|A|} {1\hspace{-3.8pt} 1}_{A}(Y)} {\frac{|B|}{|A|}\cdot \frac{1}{|B|} {1\hspace{-3.8pt} 1}_{B}(Y)}} \class{fragment}{{}= \frac{ {1\hspace{-3.8pt} 1}_{A}(Y)} { {1\hspace{-3.8pt} 1}_{B}(Y)}} \class{fragment}{{}= {1\hspace{-3.8pt} 1}_{A}(Y)} \end{align*} \]
\[ \begin{align*} \class{fragment}{{}\text{Enfin}, U \leq {1\hspace{-3.8pt} 1}_{A}(Y) \iff {1\hspace{-3.8pt} 1}_{A}(Y)=1}\class{fragment}{{}\iff Y \in A} \end{align*} \]
Dans l’exemple précédent, on a pu appliquer la méthode de rejet sans la connaissance de \(m\)
Point important: parfois la connaissance de \(m\) est difficile à obtenir, et il est préférable de ne pas l’utiliser (notamment quand les constantes de normalisation des densités sont difficiles à calculer).
Exemples: statistiques bayesiennes, modèles graphiques, etc.
Loi de Bernoulli: avec \(U_1, \ldots, U_n\) i.i.d uniformes sur \([0,1]\) (méthode d’inversion): \[ X_i \triangleq {1\hspace{-3.8pt} 1}_{\{U_i \leq p\}} \sim \mathcal{B}(p) \]
Loi binomiale: \[ \sum_{i=1}^n {1\hspace{-3.8pt} 1}_{\{U_i \leq p\}} \sim \mathcal{B}(n,p) \] en rappelant que \[ X = X_1 + \cdots + X_n \sim \mathcal{B}(n,p) \]
Note
La méthode d’inversion marche aussi, mais nécessite le calcul de l’inverse généralisée de \(F\), donc de coefficients binomiaux…
Rappel: \(\quad X \sim \mathcal{P}(\lambda) \iff \mathbb{P}(X = k) = e^{-\lambda} \dfrac{\lambda^k}{k!}\,, \quad \forall k \in \mathbb{N}.\)
Proposition 2 (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, \(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)\).
Point numérique
Méthode adaptée par numpy.random.poisson, cf. code source, qui fut proposée par D. Knuth; Source: Wikipedia
La preuve repose sur le résultat suivant:
Lemme 1 (Loi de Erlang)
Soient \(E_1, \dots, E_n\) des v.a., i.i.d. de loi exponentielle de paramètre \(\lambda >0\). La somme \(S_n=E_1+\dots+E_n\) suit une loi d’Erlang de paramètres \((n,\lambda)\), de fonction de répartition \[ F_{n,\lambda}(t) = 1 - \sum_{k=0}^{n-1} e^{-\lambda t} \frac{(\lambda t)^k}{k!}\,. \]
Preuve partielle: Transformée de Laplace de \(\mathcal{E}(\lambda)\): \(\mathbb{E}\left(e^{-tE_1}\right) = \int_0^{+\infty} e^{-t x} \lambda e^{-\lambda x} \, d x = \frac{\lambda}{\lambda+t}\)
Densité d’une loi \(\Gamma(\alpha,\beta): \quad f(x) = \tfrac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}\)
Transformée de Laplace d’une loi \(\Gamma(\alpha,\beta): \quad \mathbb{E}\left(e^{-tX}\right) = \left(\tfrac{\beta}{\beta+t}\right)^\alpha\)
Enfin, \(\mathbb{E}\left(e^{-t S_n}\right) = \left(\mathbb{E}\left(e^{-t E_1}\right)\right)^n=\left(\tfrac{\lambda}{\lambda+t}\right)^n\)
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*} \]
Du lemme précédent : \(\mathbb{P}(S_n \leq 1) = 1 - \displaystyle\sum_{k=0}^{n-1} e^{-\lambda} \dfrac{\lambda^k}{k!}\) et \(\mathbb{P}(S_{n+1} \leq 1) = 1 - \displaystyle\sum_{k=0}^{n} e^{-\lambda} \dfrac{\lambda^k}{k!}\).
Puis, \[ \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 pour \(T \triangleq \sup \{k \in \mathbb{N}^* : S_k \leq 1\}\) \[ \mathbb{P}(T=n) = \mathbb{P}(S_n \leq 1 < S_{n+1})\,. \]
Simulation