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.

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).