Lois continues

La page Fondations a posé les outils pour décrire une v.a. continue, densité, fonction de répartition, espérance et variance, indépendance, et la synthèse qui les chapeaute : une loi est un modèle paramétré $p(x;\theta )$ (§5). La page précédente en a déroulé le catalogue discret ; celle-ci fait de même côté continu, sur quatre lois usuelles, uniforme, exponentielle, normale et normale multidimensionnelle, présentées dans cet ordre, chacune un modèle $p(x;\theta )$ selon la même trame : modélisation, illustration, propriétés, puis exemple numérique. La fiche de la loi normale culmine avec le théorème central limite, le résultat qui explique pourquoi cette loi est partout. Deux sections closent ensuite le catalogue sur ce que l'on fait d'un modèle une fois choisi : en produire des tirages, par transformation d'une variable aléatoire, puis vérifier qu'il convient aux données.

Loi uniforme continue, $X\sim \mathcal{U}(a,b)$

La loi uniforme continue modélise une valeur prise « au hasard » dans un intervalle $[a,b]$. Souvent ce type de loi est choisi lorsque nous n'avons aucun a priori sur le comportement d'une variable aléatoire, à l'exception de l'intervalle des valeurs possibles. À titre d'exemple, cette loi peut servir à modéliser l'erreur de quantification d'un convertisseur analogique-numérique (l'arrondi au pas $q$ près laisse une erreur $\mathcal{U}(-q/2,q/2)$).

Modélisation

La loi uniforme continue est paramétrée par $\theta =(a,b)$, les bornes de l'intervalle. Sous la forme générique du cours, $p(x;\theta )$ désigne ici la densité $f_X(x)$.

Modèle (Loi uniforme continue)

• PDF :

$$f_X(x)=\frac{1}{b-a},x\in [a,b]$$

• CDF :

$$F_X(x)=\frac{x-a}{b-a},x\in [a,b]$$

Illustration

Figure 1.16a, PDF : rectangles de hauteur $1/(b-a)$.

Figure 1.16b, CDF : rampe linéaire sur $[a,b]$.

Propriétés

Moments

$$\begin{array}{rl}E(X)& =(a+b)/2\\ \mathrm{var}(X)& =(b-a{)}^2/12\end{array}$$

Le centre reste le même que pour la version discrète ; la dispersion croît avec la largeur de l'intervalle, au carré.

Exemple

Le code suivant montre comment tirer 100 000 erreurs de quantification au pas $q=1$, uniformes sur $[-0,5,0,5]$.

from scratch

import numpy as np
rng = np.random.default_rng(0)

u = rng.uniform(-0.5, 0.5, size=100_000)
u.mean(), u.var()   # >>> (-0.0004, 0.0832)  (≈ 0, q²/12 = 0.0833)

scipy

import scipy.stats as sta

X = sta.uniform(loc=-0.5, scale=1.0)   # loc = a, scale = b - a (et non b)
X.var()       # >>> 0.0833
X.cdf(0.25)   # >>> 0.75

Loi exponentielle, $X\sim \mathrm{Exp}(\lambda )$

La loi exponentielle modélise la durée d'attente d'un événement sans mémoire : le temps entre deux pannes, entre deux requêtes. C'est la jumelle continue de la Poisson : si les événements se comptent en Poisson de taux $\lambda$, le temps entre deux événements suit $\mathrm{Exp}(\lambda )$.

Modélisation

La loi exponentielle est paramétrée par $\theta =\lambda >0$, le taux d'occurrence ; plus il est élevé, plus les attentes sont courtes. Sa densité est la forme $p(x;\theta )$ de ce modèle.

Modèle (Loi exponentielle)

• PDF :

$$f_X(x)=\lambda e^{-\lambda x},x\ge 0$$

• CDF :

$$F_X(x)=1-e^{-\lambda x},x\ge 0$$

Illustration

Figure 1.17a, PDF : part de $\lambda$, décroît en $e^{-\lambda x}$.

Figure 1.17b, CDF : $1-e^{-\lambda x}$, croissante vers 1.

Propriétés

Moments

$$\begin{array}{rl}E(X)& =1/\lambda \\ \mathrm{var}(X)& =1/{\lambda }^2\end{array}$$

L'écart-type est égal à la moyenne : la loi est fortement asymétrique, elle servira de loi « pas du tout Gaussienne » pour mettre le TCL à l'épreuve en fin de fiche normale.

Preuve

Par intégration par parties, $E(X)={\int }_0^{\mathrm{\infty }}x\lambda e^{-\lambda x}\mathrm{d}x=1/\lambda$, et $E(X^2)=2/{\lambda }^2$, d'où $\mathrm{var}(X)=2/{\lambda }^2-1/{\lambda }^2=1/{\lambda }^2$ par König-Huygens. ∎

L'absence de mémoire

Proposition, Absence de mémoire (exponentielle)

Si $X\sim \mathrm{Exp}(\lambda )$, alors pour tous $s,t\ge 0$ :

$$P(X>s+t\mid X>s)=P(X>t).$$

Une machine qui fonctionne depuis $s$ heures n'est ni plus ni moins « due » pour une panne qu'une machine neuve. Il peut être démontré que l'exponentielle est la seule loi continue à posséder cette propriété ; son analogue discret est la géométrique.

Démonstration

$P(X>s+t\mid X>s)=P(X>s+t)/P(X>s)=e^{-\lambda (s+t)}/e^{-\lambda s}=e^{-\lambda t}=P(X>t)$. ∎

Exemple

Le code suivant montre comment tirer 100 000 durées entre pannes d'un serveur de taux $\lambda =0,5$ panne/h, puis vérifier l'absence de mémoire.

from scratch

import numpy as np
rng = np.random.default_rng(0)

t = rng.exponential(scale=2.0, size=100_000)   # scale = 1/λ
t.mean(), t.var()   # >>> (1.994, 4.016)  (≈ 2, 4)
np.mean(t > 3)      # >>> 0.2226  (≈ e^{-1.5})

scipy

import scipy.stats as sta

X = sta.expon(scale=2.0)
X.sf(3)             # >>> 0.2231
X.sf(5) / X.sf(2)   # >>> 0.2231  (= P(X > 3) : sans mémoire)

Loi normale, $X\sim \mathcal{N}(\mu ,{\sigma }^2)$

La loi normale modélise le bruit de mesure, et plus généralement toute grandeur résultant de l'accumulation de petits effets indépendants. Sa densité est la célèbre courbe en cloche, symétrique.

Modélisation

La loi normale est paramétrée par $\theta =(\mu ,{\sigma }^2)$ : $\mu$ positionne le centre de la cloche, ${\sigma }^2$ règle sa largeur ; sa densité est la forme $p(x;\theta )$ de ce modèle. Convention absolue du cours : le second argument de $\mathcal{N}(\mu ,{\sigma }^2)$ est la variance, pas l'écart-type.

Modèle (Loi normale)

• PDF :

$$f_X(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{(x-\mu {)}^2}{2{\sigma }^2})$$

• CDF :

$$F_X(x)=\mathrm{\Phi }\left(\frac{x-\mu }{\sigma }\right)$$

où $\mathrm{\Phi }$ est la fonction de répartition de la normale centrée réduite (ci-dessous), pas de forme fermée élémentaire : tout calcul se ramène à la table (ou à scipy) de $\mathrm{\Phi }$, via un changement d'échelle.

Illustration

Figure 1.18a, PDF : $\mu$ déplace, ${\sigma }^2$ élargit.

Figure 1.18b, CDF : une rampe en S, centrée en $\mu$.

Propriétés

Moments

$$\begin{array}{rl}E(X)& =\mu \\ \mathrm{var}(X)& ={\sigma }^2\end{array}$$

Remarquons que les paramètres de la loi sont exactement ses deux premiers moments.

La standardisation

Toutes les Gaussiennes se déduisent d'une seule, $\mathcal{N}(0,1)$, par un simple changement d'échelle : soustraire $\mu$, diviser par $\sigma$. Tout calcul de probabilité Gaussienne se ramène ainsi à une table (ou une fonction) unique, $\mathrm{\Phi }$.

Proposition, Standardisation (centrage-réduction)

Si $X\sim \mathcal{N}(\mu ,{\sigma }^2)$, alors

$$Z=\frac{X-\mu }{\sigma }\sim \mathcal{N}(0,1),$$

la loi normale centrée réduite, de fonction de répartition notée $\mathrm{\Phi }$. Tout calcul Gaussien s'y ramène :

$$P(X\le x)=\mathrm{\Phi }\left(\frac{x-\mu }{\sigma }\right).$$

La règle 68–95–99,7

La standardisation fige quelques probabilités de référence, valables pour n'importe quelle Gaussienne.

Écart au centre	Intervalle	$P(\mu -k\sigma \le X\le \mu +k\sigma )$
$k=1$	$[\mu -\sigma ,\mu +\sigma ]$	$68,3\mathrm{\%}$
$k=1,96$	$[\mu -1,96\sigma ,\mu +1,96\sigma ]$	$95,0\mathrm{\%}$
$k=2$	$[\mu -2\sigma ,\mu +2\sigma ]$	$95,4\mathrm{\%}$
$k=3$	$[\mu -3\sigma ,\mu +3\sigma ]$	$99,7\mathrm{\%}$

C'est la règle mnémotechnique 68–95–99,7. Le quantile $1,96$ est celui des 95 % exacts, celui des intervalles de confiance usuels.

Figure 1.19, La règle 68–95–99,7 : quelle que soit la Gaussienne, environ 68 % de la probabilité vit à $\pm \sigma$ du centre, 95 % à $\pm 2\sigma$, 99,7 % à $\pm 3\sigma$. Une observation à plus de $3\sigma$ est un événement à 3 pour 1 000, le cœur du raisonnement des tests du chapitre 3 (script figures/ch1/03-normale-68-95-997.py).

Le théorème central limite

Reste la question soulevée en tête de fiche : pourquoi cette loi serait-elle si répandue, alors que le réel n'est pas toujours en cloche (le marathon) ? La réponse tient dans le théorème le plus célèbre des probabilités. Il porte sur la somme de v.a. i.i.d., abréviation d'indépendantes et identiquement distribuées : les $X_i$ suivent toutes la même loi, et sont mutuellement indépendantes. C'est le modèle par défaut de $n$ mesures répétées d'un même phénomène, sans influence des unes sur les autres, et l'hypothèse de travail de tout le cours.

Deux façons pour une suite de variables aléatoires de « se rapprocher » d'une limite

Une suite de v.a. $Y_n$ converge en probabilité vers une constante $c$ si, pour tout $\epsilon >0$, $P(|Y_n-c|>\epsilon )\to 0$ quand $n\to \mathrm{\infty }$ : la loi de $Y_n$ se concentre sur $c$, qui doit être un nombre fixe. Elle converge en loi vers une v.a. $Y$ si sa fonction de répartition converge vers celle de $Y$ en tout point de continuité ; $Y$ peut ici être elle-même aléatoire, avec une loi non dégénérée, c'est ce mode de convergence qu'utilise le théorème ci-dessous, la limite étant une loi normale, pas un nombre.

Théorème central limite (TCL)

Soient $X_1,\dots ,X_n$ des v.a. i.i.d. d'espérance $\mu$ et de variance ${\sigma }^2$ finie, quelle que soit leur loi, et $S_n=X_1+\cdots +X_n$ leur somme. Alors, quand $n\to \mathrm{\infty }$, la somme standardisée converge en loi vers la Gaussienne centrée réduite :

$$\frac{S_n-n\mu }{\sigma \sqrt{n}}\underset{n\to \mathrm{\infty }}{\overset{\text{loi}}{\to }}\mathcal{N}(0,1),$$

ce qui s'utilise en pratique sous la forme : pour $n$ assez grand,

$$S_n\approx \mathcal{N}(n\mu ,n{\sigma }^2).$$

L'énoncé est admis. Son point remarquable est l'universalité : quelle que soit la loi de départ, dé, exponentielle, Bernoulli, sommer suffit à faire émerger la Gaussienne. D'où la justification de la phrase d'ouverture de cette fiche : toute grandeur qui accumule un grand nombre de petits effets indépendants (un bruit de capteur, une somme d'erreurs élémentaires) est presque une somme de v.a. i.i.d., donc approximativement Gaussienne par TCL. C'est cet argument qui motive, dans tout le cours, la modélisation Gaussienne des composantes de bruit.

Exemple de convergence : sommer des exponentielles

Le théorème se met à l'épreuve sur la loi la plus asymétrique du catalogue, $\mathrm{Exp}(1)$ ($\mu =1$, ${\sigma }^2=1$) : rien ne ressemble moins à une cloche que sa densité.

Figure 1.20, Gaussianisation de la somme $S_N$ pour une loi de départ $\mathrm{Exp}(1)$, très asymétrique : à $N=1$ la prédiction Gaussienne est franchement fausse (et doit l'être), à $N=10$ elle est déjà bonne, à $N=100$ elle est indiscernable de l'histogramme (script figures/ch1/03-tcl-convergence.py).

Conséquence : la loi de la moyenne

L'analyse de données ne somme pas ses mesures, elle les moyenne. Le théorème s'y transpose sans calcul supplémentaire : la moyenne empirique ${\bar{X}}_n=S_n/n$ n'est que la somme ramenée à l'échelle d'une observation, et diviser par $n$ ne change pas la forme d'une loi.

Corollaire, loi de la moyenne empirique

Sous les hypothèses du théorème central limite, pour $n$ assez grand,

$${\bar{X}}_n=\frac{1}{n}\sum _{i=1}^n{X}_i\approx \mathcal{N}(\mu ,\frac{{\sigma }^2}{n}).$$

Les deux paramètres se lisent directement sur ceux de $S_n$ : $E({\bar{X}}_n)=n\mu /n=\mu$, et $\mathrm{var}({\bar{X}}_n)=n{\sigma }^2/n^2={\sigma }^2/n$ par la règle $\mathrm{var}(aX)=a^2\mathrm{var}(X)$. Le contraste avec la somme est instructif. La dispersion de la somme croît en $\sigma \sqrt{n}$ : additionner des aléas les accumule. Celle de la moyenne fond en $\sigma /\sqrt{n}$ : la moyenne vise juste en espérance et se resserre autour de sa cible à mesure que $n$ grandit. Diviser par 10 l'incertitude d'une moyenne exige donc 100 fois plus de mesures, une contrainte que tout ingénieur finit par rencontrer.

C'est sous cette forme que le TCL servira dans tout le cours : la Gaussienne est la loi des moyennes. Le chapitre 2 en fait son point de départ, la moyenne empirique y devenant le premier estimateur rencontré.

Loi normale multidimensionnelle, $\mathbf{X}\sim \mathcal{N}(\mathit{\mu },\mathbf{C})$

Une mesure isolée n'est pas toujours la situation la plus riche : deux capteurs relevés ensemble, deux coordonnées d'une position, ou plus généralement un vecteur de $n$ mesures forment un objet aléatoire à part entière. La Gaussienne se généralise à un tel vecteur $\mathbf{X}=(X_1,\dots ,X_n)$.

Modélisation

La loi normale multidimensionnelle est paramétrée par $\theta =(\mathit{\mu },\mathbf{C})$ : $\mathit{\mu }\in {\mathbb{R}}^n$ le vecteur des espérances, et $\mathbf{C}$, une matrice $n\times n$ symétrique définie positive, la matrice de covariance ($C_{ii}$ la variance de $X_i$, $C_{ij}$ la covariance entre $X_i$ et $X_j$). Sa densité est la forme $p(\mathbf{x};\theta )$ de ce modèle. Elle n'admet pas de CDF sous forme fermée.

Modèle (Loi normale multidimensionnelle)

• PDF :

$$f_{\mathbf{X}}(\mathbf{x})=\frac{1}{(2\pi {)}^{n/2}|\mathbf{C}{|}^{1/2}}\exp (-\frac{1}{2}(\mathbf{x}-\mathit{\mu }{)}^T{\mathbf{C}}^{-1}(\mathbf{x}-\mathit{\mu }))$$

Illustration

Pour $n=2$, avec des variances ${\sigma }_1^2,{\sigma }_2^2$ et une corrélation $\rho$,

$$\mathbf{C}=\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right).$$

Ce que $\mathbf{C}$ change à la forme du nuage se lit directement sur deux simulations en dimension 2 :

Figure 1.21a, $X_1$ et $X_2$ non corrélées, même dispersion sur chaque axe (${\sigma }_1={\sigma }_2=1$), 500 tirages simulés : le nuage est rond.

Figure 1.21b, $X_1$ et $X_2$ corrélées ($\rho =0,8$), même dispersion sur chaque axe : le nuage s'incline et s'étire selon la diagonale (script figures/ch1/03-normale-2d.py).

Le nuage rond correspond à $\rho =0$, soit une matrice diagonale $\mathbf{C}=\mathbf{I}$ ; le nuage incliné correspond à $\rho =0,8$, dont les termes hors diagonale non nuls inclinent l'ellipse selon la diagonale.

Propriétés

Moments

$$\begin{array}{rl}E(\mathbf{X})& =\mathit{\mu }\\ \mathrm{cov}(\mathbf{X})& =\mathbf{C}\end{array}$$

Le cas non corrélé

Si $\mathbf{C}$ est diagonale, $\mathbf{C}=\mathrm{diag}({\sigma }_1^2,\dots ,{\sigma }_n^2)$, alors les composantes $X_1,\dots ,X_n$ sont indépendantes, et la densité se factorise :

$$f_{\mathbf{X}}(\mathbf{x})=\prod _{i=1}^n{f}_{X_i}(x_i),X_i\sim \mathcal{N}({\mu }_i,{\sigma }_i^2).$$

Pour une Gaussienne, absence de corrélation et indépendance coïncident, un résultat qui ne vaut pas pour une loi quelconque. Un bruit blanc ($\mathbf{C}={\sigma }^2\mathbf{I}$) n'est alors rien d'autre que $n$ bruits Gaussiens scalaires indépendants, empilés dans un vecteur.

Exemple

Le code suivant montre comment évaluer la densité de deux capteurs bruités corrélés à $\rho =0,8$, tous deux d'écart-type unité, puis en tirer 100 000 réalisations.

from scratch

import numpy as np

def densite_normale_2d(x, mu, C):
    """Densité de N(mu, C) en dimension 2."""
    ecart = x - mu
    inv_C = np.linalg.inv(C)
    det_C = np.linalg.det(C)
    exposant = -0.5 * ecart @ inv_C @ ecart
    return np.exp(exposant) / (2 * np.pi * np.sqrt(det_C))

mu = np.array([0.0, 0.0])
C = np.array([[1.0, 0.8], [0.8, 1.0]])
densite_normale_2d(np.array([0.0, 0.0]), mu, C)   # >>> 0.2653  (le pic)

scipy

import scipy.stats as sta

X = sta.multivariate_normal(mean=[0.0, 0.0], cov=[[1.0, 0.8], [0.8, 1.0]])
X.pdf([0.0, 0.0])   # >>> 0.2653

ech = X.rvs(size=100_000, random_state=0)
np.cov(ech.T)       # >>> [[1.00, 0.80], [0.80, 1.01]]  (≈ C)

Au chapitre 2 : ce modèle, plusieurs mesures bruitées et corrélées à la fois, sera au cœur de la dérivation des moindres carrés.

Récapitulatif : une loi = un modèle paramétré

Le catalogue est complet. Reprenons la lecture posée en fondations, §5, une loi est un modèle paramétré $p(x;\theta )$, et appliquons-la aux huit lois des deux catalogues :

Loi	$\theta$	$E(X)$	$\mathrm{var}(X)$	Situation-type
Bernoulli $\mathcal{B}(p)$	$p$	$p$	$p(1-p)$	conversion d'un test A/B
Binomiale $\mathcal{B}(n,p)$	$(n,p)$	$np$	$np(1-p)$	bits erronés sur une trame
Géométrique $\mathcal{G}(p)$	$p$	$1/p$	$(1-p)/p^2$	attente du premier succès
Poisson $\mathcal{P}(\lambda )$	$\lambda$	$\lambda$	$\lambda$	requêtes par seconde
Uniforme discrète $\mathcal{U}\{a,b\}$	$(a,b)$	$(a+b)/2$	$(n^2-1)/12$	le dé
Uniforme $\mathcal{U}(a,b)$	$(a,b)$	$(a+b)/2$	$(b-a{)}^2/12$	erreur de quantification
Exponentielle $\mathrm{Exp}(\lambda )$	$\lambda$	$1/\lambda$	$1/{\lambda }^2$	temps entre deux pannes
Normale $\mathcal{N}(\mu ,{\sigma }^2)$	$(\mu ,{\sigma }^2)$	$\mu$	${\sigma }^2$	bruit de mesure

Huit situations, une seule idée : chaque ligne est un modèle $p(x;\theta )$, et modéliser, c'est choisir la ligne. Chacune de ces lignes s'interroge avec l'interface unique de scipy.stats (fondations, §5.3) : le même code fonctionne pour les huit lois en changeant la première ligne, sta.bernoulli(p), sta.binom(n,p), sta.geom(p), sta.poisson(lam), sta.randint(a, b+1), sta.uniform(loc=a, scale=b-a), sta.expon(scale=1/lam), sta.norm(loc=mu, scale=sigma).

Transformation d'une variable aléatoire

Les sections précédentes décrivent une variable aléatoire dont la loi est donnée. Deux besoins inverses se présentent en pratique : déduire la loi d'une quantité calculée à partir d'une autre, $Y=g(X)$, et produire des tirages d'une loi voulue à partir du seul générateur dont dispose une machine, un tirage uniforme sur $[0,1]$. Les deux relèvent de la même opération, transformer une variable aléatoire.

Changement de variable

Proposition, Changement de variable

Soit $X$ une v.a. continue de densité $f_X$ et $g$ une fonction strictement monotone et dérivable. La densité de $Y=g(X)$ est

$$f_Y(y)=f_X(g^{-1}(y))|\frac{\mathrm{d}{g}^{-1}}{\mathrm{d}y}(y)|,$$

où

• $g^{-1}$ est la fonction réciproque de $g$,
• le facteur $|\mathrm{d}{g}^{-1}/\mathrm{d}y|$ corrige la dilatation locale de l'axe.

Démonstration

Pour $g$ croissante, l'événement $\{Y\le y\}$ est l'événement $\{X\le g^{-1}(y)\}$, donc les fonctions de répartition se répondent :

$$F_Y(y)=F_X(g^{-1}(y)).$$

Dériver les deux membres par rapport à $y$ donne la densité annoncée. Pour $g$ décroissante, l'inégalité se retourne et fait apparaître un signe moins, que la valeur absolue absorbe. $◼$

Le cas de loin le plus fréquent est la transformation affine, $Y=aX+b$, pour laquelle $g^{-1}(y)=(y-b)/a$ et le facteur correctif vaut $1/|a|$. Ses moments se lisent sans calcul, par la linéarité de l'espérance et la règle sur la variance : $E(Y)=aE(X)+b$ et $\mathrm{var}(Y)=a^2\mathrm{var}(X)$. Appliquée à la Gaussienne, elle produit la standardisation, qui ramène toute cloche à la cloche centrée réduite.

Méthode de la transformation inverse

Un ordinateur ne sait engendrer qu'une seule loi, l'uniforme sur $[0,1]$. Toutes les autres s'en déduisent, et la fonction de répartition fournit le chemin.

Théorème, Transformation inverse

Soit $F_X$ une fonction de répartition continue et strictement croissante, et $U\sim \mathcal{U}(0,1)$. La variable aléatoire

$$X=F_X^{-1}(U)$$

admet $F_X$ pour fonction de répartition.

Démonstration

$F_X$ étant croissante, l'événement $\{F_X^{-1}(U)\le x\}$ est l'événement $\{U\le F_X(x)\}$. Or $U$ est uniforme sur $[0,1]$, donc $P(U\le u)=u$ pour tout $u$ de cet intervalle. En prenant $u=F_X(x)$ :

$$P(X\le x)=P(U\le F_X(x))=F_X(x).◼$$

Le mécanisme se comprend mieux en le regardant opérer. Tirer $u$ uniformément, c'est tirer uniformément une probabilité cumulée ; la renvoyer sur l'axe des $x$ par $F_X^{-1}$ concentre les valeurs là où la fonction de répartition monte vite, c'est-à-dire là où la densité est forte.

Illustration

Si le théorème est exact, des $u$ régulièrement répartis sur l'axe vertical doivent ressortir resserrés là où la densité visée est grande, et un grand nombre de tirages doit reconstituer cette densité.

Figure 1.22a, La transformation inverse en action sur la loi exponentielle ($\lambda =1$) : cinq valeurs de $u$ renvoyées sur l'axe des $x$ par $F_X^{-1}$ (script figures/ch1/03-transformation-inverse.py).

Figure 1.22b, Dix mille tirages produits par cette méthode, confrontés à la densité visée.

À gauche, les cinq $u$ sont régulièrement espacés sur l'axe vertical, et les cinq $x$ qu'ils produisent ne le sont pas : quatre d'entre eux tombent en deçà de $1,5$, parce que c'est là que la fonction de répartition gagne le plus de hauteur. La méthode transporte donc l'uniformité des probabilités cumulées vers la non-uniformité des valeurs. À droite, l'histogramme de dix mille tirages épouse la densité exponentielle, et leur moyenne vaut $1,02$ pour une espérance théorique de $1/\lambda =1$.

Le cas discret

Une fonction de répartition discrète progresse par paliers et n'est donc pas inversible au sens strict. La construction se transpose en retenant la première valeur dont la probabilité cumulée dépasse $u$ :

$$X=min\{x_k:F_X(x_k)\ge U\}.$$

Géométriquement, l'axe $[0,1]$ est découpé en segments dont les longueurs sont les probabilités $P(X=x_k)$ ; le tirage $u$ tombe dans l'un d'eux, et désigne la valeur correspondante. Une valeur deux fois plus probable occupe un segment deux fois plus long, donc est retenue deux fois plus souvent.

Implémentation

La méthode s'écrit en trois étapes, identiques dans les deux cas :

1. tirer $u$ uniformément sur $[0,1]$, autant de fois que de valeurs souhaitées ;
2. inverser la fonction de répartition, analytiquement si elle s'inverse, par recherche dans les probabilités cumulées sinon ;
3. vérifier l'échantillon obtenu contre la loi visée, moyenne et variance empiriques, ou histogramme.

from scratch

import numpy as np
rng = np.random.default_rng(2026)

M, lam = 10_000, 1.0

# continuous case: F(x) = 1 - exp(-lam x) inverts in closed form
u = rng.random(M)
x = -np.log(1.0 - u) / lam
x.mean()                                  # >>> 1.0202  (expected: 1 / lam = 1)

# discrete case: a loaded die, the first value whose cumulated mass exceeds u
pmf = np.array([0.10, 0.10, 0.10, 0.10, 0.10, 0.50])
cdf = np.cumsum(pmf)
faces = np.searchsorted(cdf, rng.random(M)) + 1
np.mean(faces == 6)                       # >>> 0.4945  (expected: 0.50)

scipy

import scipy.stats as sta

X = sta.expon(scale=1 / lam)
X.rvs(size=M, random_state=2026).mean()   # >>> 1.0047  (expected: 1 / lam = 1)

X.ppf(0.5)                                # >>> 0.6931  (the median, ln 2 / lam)

La méthode rvs de scipy.stats (section 5.3 des fondations) ne fait rien d'autre : elle applique $F_X^{-1}$, exposée par ailleurs sous le nom ppf (percent point function), à des tirages uniformes.

Limites

La méthode exige de savoir inverser $F_X$, ce qui n'est pas toujours possible. La loi normale en est l'exemple central : sa fonction de répartition $\mathrm{\Phi }$ n'a pas de forme fermée élémentaire, et son inverse non plus. Deux voies contournent l'obstacle. La première inverse $\mathrm{\Phi }$ numériquement, c'est ce que fait sta.norm.ppf. La seconde emprunte une transformation différente, la transformation de Box-Muller, qui produit deux Gaussiennes centrées réduites indépendantes à partir de deux uniformes indépendantes :

$$Z_1=\sqrt{-2\ln U_1}\cos (2\pi U_2),$$$$Z_2=\sqrt{-2\ln U_1}\sin (2\pi U_2).$$

Sa démonstration relève du changement de variable en dimension deux ; sa mise en œuvre fait l'objet du projet Les lois usuelles en NumPy, qui recode ainsi les tirages de tout le catalogue du chapitre.

Adéquation d'un modèle

Le catalogue fournit huit modèles ; rien n'a encore dit comment vérifier que celui retenu convient aux données. Comment confronter une loi $p(x;\theta )$ aux mesures qu'elle prétend décrire ? Deux lectures graphiques suffisent à trancher la plupart des cas, avant tout calcul.

Définition

Définition, Diagramme quantile-quantile

Le diagramme quantile-quantile (QQ-plot) porte les quantiles du modèle en abscisse et les quantiles observés en ordonnée. Pour $N$ mesures triées $x_{(1)}\le \cdots \le x_{(N)}$, il trace les $N$ points

$$(F^{-1}(\frac{i-0,5}{N};\theta ),x_{(i)}),i=1,\dots ,N,$$

où

• $F^{-1}(\cdot ;\theta )$ est la fonction quantile du modèle, l'inverse de sa fonction de répartition,
• $x_{(i)}$ est la $i$-ème plus petite mesure,
• $(i-0,5)/N$ est la fraction de l'échantillon située sous cette mesure.

Si le modèle décrit les données, les points s'alignent sur la première bissectrice.

La lecture est immédiate : un point au-dessus de la bissectrice signale une mesure plus grande que ce que le modèle prévoit à ce rang, un point en dessous une mesure plus petite. Un décrochage aux deux extrémités, en arc, révèle une asymétrie que le modèle ne reproduit pas.

Illustration

L'histogramme des chronos de marathon de la figure 1.4 s'étirait vers la droite. Si la loi normale, symétrique, décrit mal cette asymétrie, l'écart doit rester discret sur l'histogramme et éclater sur le diagramme quantile-quantile.

Figure 1.23a, Chronos de $1998$ marathoniens simulés et loi normale de même moyenne et de même écart-type qu'eux, $4,08$ h et $0,68$ h (script figures/ch1/03-adequation-marathon.py).

Figure 1.23b, Diagramme quantile-quantile des mêmes chronos contre cette loi normale ; la bissectrice correspond à un modèle exact.

À gauche, la Gaussienne suit correctement le corps de la distribution, et l'œil hésite à la rejeter. À droite, les points quittent la bissectrice aux deux extrémités et dessinent un arc, la signature d'une loi asymétrique confrontée à un modèle symétrique. Le désaccord se chiffre : le modèle prévoit $5$ coureurs au-delà de $6$ h parmi les $1998$, il y en a $28$ ; au-delà de $6$ h $30$, il en prévoit moins d'un et il y en a $12$. L'histogramme dissimule ces écarts parce qu'ils portent sur de faibles effectifs ; le diagramme quantile-quantile les expose parce qu'il compare les données au modèle rang par rang.

Confusion classique

« La densité superposée épouse l'histogramme, donc le modèle convient. » Non : l'histogramme est dominé par le corps de la distribution, là où se concentrent les effectifs. Les queues, qui portent les événements rares et souvent les plus coûteux, y sont pratiquement invisibles.

Implémentation

Le diagramme se construit en trois lignes, et scipy en fournit une version toute faite :

from scratch

import numpy as np
import scipy.stats as sta
rng = np.random.default_rng(2026)

# the marathon times of figure 1.4, simulated from a shifted lognormal
chronos = 2.0 + rng.lognormal(mean=np.log(2.0), sigma=0.32, size=2_000)
chronos = chronos[chronos < 7.5]

x = np.sort(chronos)                              # the N measurements, sorted
p = (np.arange(1, x.size + 1) - 0.5) / x.size     # plotting positions
q = sta.norm.ppf(p, x.mean(), x.std())            # quantiles of the model
# then plot x against q, together with the first bisector

scipy

import matplotlib.pyplot as plt
import scipy.stats as sta

sta.probplot(chronos, dist="norm", plot=plt)
# quantiles are standardised here: the reference is the least-squares line,
# not the first bisector

Ces deux lectures restent qualitatives : elles montrent un désaccord sans le quantifier, et deux personnes peuvent lire le même diagramme différemment. Le chapitre 3 transforme ce jugement en décision chiffrée, en traitant l'adéquation d'un modèle comme un problème de détection.