Étudiez plusieurs variables à la fois

Une station météo rennaise enregistre chaque week-end le nombre de jours de pluie et la présence de vent fort. Prises séparément, ces deux variables sont faciles à décrire. Mais la question intéressante est ailleurs : quand il y a du vent, pleut-il davantage ? Pour y répondre, il faut étudier les deux variables ensemble.

Ce chapitre généralise tout ce que vous avez vu à plusieurs variables : loi conjointe, lois marginales, indépendance, covariance, et enfin les vecteurs gaussiens, omniprésents en statistique et en apprentissage automatique. Le fil rouge sera cette station météo, puis deux capteurs de température.

Loi conjointe et lois marginales

Un vecteur aléatoire $\mathbf{X} = (X_1, \dots, X_n)$ est une fonction de $\Omega$ dans $\mathbb{R}^n$ : $n$ variables aléatoires définies sur la même expérience. Pour deux variables, on parle de couple $(X, Y)$.

Le cas discret : un tableau

La loi conjointe d’un couple discret est le tableau des probabilités $P(X = x_i, Y = y_j)$ (la virgule se lit « et »). Pour notre station météo, $X$ est le nombre de jours de pluie sur le week-end (0, 1 ou 2) et $Y$ vaut 1 s’il y a du vent fort, 0 sinon. Supposons que l’historique donne :

	$Y = 0$ (pas de vent fort)	$Y = 1$ (vent fort)	Loi de $X$
$X = 0$	0,30	0,05	0,35
$X = 1$	0,20	0,10	0,30
$X = 2$	0,15	0,20	0,35
Loi de $Y$	0,65	0,35	1

Les lois marginales s’obtiennent en sommant les lignes ou les colonnes, d’où leur nom (elles s’écrivent dans les marges) :

$$P(X = x_i) = \sum_j P(X = x_i, Y = y_j)$$

C’est simplement la formule des probabilités totales du chapitre 2, avec la partition des valeurs de $Y$.

La loi conjointe au centre ; chaque loi marginale est la somme d'une ligne ou d'une colonne, écrite « dans la marge ».

La loi conditionnelle de $X$ sachant $Y = y_j$ s’obtient en normalisant une colonne : $P(X = x_i \mid Y = y_j) = P(X = x_i, Y = y_j) / P(Y = y_j)$. Les week-ends venteux, $P(X = 2 \mid Y = 1) = 0{,}20 / 0{,}35 \approx 0{,}57$, contre $0{,}35$ tous week-ends confondus : le vent fort s’accompagne plus souvent de pluie.

Conditionner, c'est ne garder qu'une colonne du tableau et la renormaliser : les week-ends venteux, deux jours de pluie deviennent le cas le plus probable (0,57 contre 0,35).

Le cas continu : une densité conjointe

Pour un couple continu, la loi conjointe est une densité $f(x, y) \geq 0$ d’intégrale 1 sur $\mathbb{R}^2$, et $P((X, Y) \in D) = \iint_D f(x, y) \, dx \, dy$ est un volume sous la surface. Les formules discrètes se transposent en remplaçant les sommes par des intégrales :

$$f_X(x) = \int_{-\infty}^{+\infty} f(x, y) \, dy, \qquad f_{X \mid Y = y}(x) = \frac{f(x, y)}{f_Y(y)}$$

La première donne la densité marginale de $X$ ; la seconde, la densité conditionnelle de $X$ sachant $Y = y$ (pour $f_Y(y) \gt 0$).

Une densité conjointe vue de dessus (plus c'est foncé, plus c'est probable). Les marginales sont ses « ombres » sur chaque axe ; une coupe verticale, renormalisée, donne la densité conditionnelle de Y sachant X = x.

À retenir

La loi conjointe contient toute l’information sur le couple. Les lois marginales s’en déduisent en sommant (ou en intégrant) sur l’autre variable ; l’inverse est faux : les marginales seules ne permettent pas de reconstruire la conjointe.

Indépendance de variables aléatoires

Deux variables $X$ et $Y$ sont indépendantes si tout événement portant sur $X$ est indépendant de tout événement portant sur $Y$. En pratique, on vérifie que la loi conjointe est le produit des marginales :

• cas discret : $P(X = x_i, Y = y_j) = P(X = x_i) \, P(Y = y_j)$ pour tous $i$, $j$ ;
• cas continu : $f(x, y) = f_X(x) \, f_Y(y)$ pour tous $x$, $y$ (ou de façon équivalente, avec les fonctions de répartition).

Pour la station météo, $P(X = 0, Y = 1) = 0{,}05$, alors que $P(X = 0) \, P(Y = 1) = 0{,}35 \times 0{,}35 \approx 0{,}12$ : pluie et vent ne sont pas indépendants. Une seule case qui ne factorise pas suffit.

Somme de variables indépendantes

Si $X$ et $Y$ sont indépendantes, la loi de $Z = X + Y$ se calcule par un produit de convolution. Dans le cas continu :

$$f_Z(z) = \int_{-\infty}^{+\infty} f_X(z - y) \, f_Y(y) \, dy$$

Pour que la somme vaille $z$, il faut que $Y$ prenne une valeur $y$ et $X$ la valeur $z - y$ ; on intègre sur toutes les valeurs de $y$. Dans le cas discret, l’intégrale devient une somme : c’est exactement ce que vous avez fait au chapitre 1 en comptant les couples de dés de somme 7. La convolution est aussi l’opération de base du filtrage en traitement du signal et des réseaux de neurones convolutifs.

La somme de deux variables indépendantes « étale » les deux lois : la durée totale U + T est plus large que le trajet seul et n'est plus exactement normale.

Covariance et corrélation

La loi conjointe est parfois lourde. Comme l’espérance et la variance résument une variable, la covariance résume la façon dont deux variables varient ensemble :

$$\mathrm{Cov}(X, Y) = E\big[(X - E[X])(Y - E[Y])\big] = E[XY] - E[X] \, E[Y]$$

Elle est positive si $X$ et $Y$ ont tendance à être grandes (ou petites) en même temps, négative si l’une a tendance à être grande quand l’autre est petite. On a $\mathrm{Cov}(X, X) = \mathrm{Var}(X)$. Son défaut : elle dépend des unités. On la normalise donc en coefficient de corrélation :

$$\rho(X, Y) = \frac{\mathrm{Cov}(X, Y)}{\sigma_X \, \sigma_Y} \in [-1\,;\,1]$$

$\rho = \pm 1$ signifie que $Y$ est une fonction affine de $X$ ; $\rho$ mesure donc l’intensité d’une relation linéaire.

La covariance est la moyenne du produit des écarts aux moyennes. Quand la plupart des points tombent dans les quadrants « + », elle est positive.

ρ mesure l'alignement autour d'une droite. Les deux derniers nuages ont une corrélation nulle alors que X et Y sont clairement liés : c'est le piège ci-dessous.

Exemple

Pour la station météo : $E[X] = 0 \times 0{,}35 + 1 \times 0{,}30 + 2 \times 0{,}35 = 1$ et $E[Y] = 0{,}35$. Le produit $XY$ n’est non nul que si $Y = 1$, donc $E[XY] = 1 \times 0{,}10 + 2 \times 0{,}20 = 0{,}5$. D’où :

$$\mathrm{Cov}(X, Y) = 0{,}5 - 1 \times 0{,}35 = 0{,}15$$

Avec $\mathrm{Var}(X) = E[X^2] - 1 = (0{,}30 + 4 \times 0{,}35) - 1 = 0{,}7$ et $\mathrm{Var}(Y) = 0{,}35 \times 0{,}65 = 0{,}2275$, on trouve $\rho = 0{,}15 / \sqrt{0{,}7 \times 0{,}2275} \approx 0{,}38$ : une association positive, modérée.

La covariance intervient dans la variance d’une somme :

$$\mathrm{Var}(X + Y) = \mathrm{Var}(X) + \mathrm{Var}(Y) + 2\,\mathrm{Cov}(X, Y)$$

Si $X$ et $Y$ sont indépendantes, alors $E[XY] = E[X] \, E[Y]$, donc $\mathrm{Cov}(X, Y) = 0$ et les variances s’additionnent. Revenons au trajet du chapitre 5 : l’attente du bus $U \sim \mathcal{U}[0, 10]$ et le trajet $T \sim \mathcal{N}(30, 4^2)$ sont indépendants, donc la durée totale a pour espérance $5 + 30 = 35$ minutes et pour variance $100/12 + 16 \approx 24{,}3$, soit un écart-type de 4,9 minutes. Sa loi exacte (une convolution, qui n’est pas normale) donne un quantile à 95 % d’environ 43,1 minutes : partez 44 minutes avant votre cours.

Piège

Indépendance implique covariance nulle, mais la réciproque est fausse. Prenez $X$ uniforme sur $[-1\,;\,1]$ et $Y = X^2$ : $Y$ est entièrement déterminée par $X$, et pourtant $\mathrm{Cov}(X, Y) = E[X^3] - E[X] \, E[X^2] = 0 - 0 = 0$ par symétrie. La covariance ne détecte que les liens linéaires.

Un autre contre-exemple classique : un point $(X, Y)$ tiré uniformément dans le disque unité. Par symétrie, $\mathrm{Cov}(X, Y) = 0$. Mais si $X = 0{,}9$, alors $Y$ est forcément entre $-0{,}44$ et $0{,}44$ : connaître $X$ renseigne sur $Y$. Formellement, $f(0{,}9\,;\,0{,}9) = 0$ (le point est hors du disque) alors que $f_X(0{,}9) \, f_Y(0{,}9) \gt 0$. Les variables sont non corrélées mais dépendantes.

Vecteur moyen et matrice de covariance

Avec $n$ variables, on range les espérances dans le vecteur moyen $\boldsymbol{\mu} = E[\mathbf{X}] = (E[X_1], \dots, E[X_n])$, et toutes les covariances dans la matrice de covariance :

$$\Sigma = E\big[(\mathbf{X} - \boldsymbol{\mu})(\mathbf{X} - \boldsymbol{\mu})^T\big], \qquad \Sigma_{ij} = \mathrm{Cov}(X_i, X_j)$$

Sa diagonale contient les variances. Elle est symétrique et semi-définie positive (pour tout vecteur $\mathbf{a}$, $\mathbf{a}^T \Sigma \, \mathbf{a} \geq 0$, car c’est la variance de $\mathbf{a}^T \mathbf{X}$). Si les composantes sont indépendantes, elle est diagonale. Son inverse, quand il existe, s’appelle la matrice de précision.

Tout son intérêt apparaît avec les transformations affines. Si $\mathbf{Y} = A\mathbf{X} + \mathbf{b}$, avec $A$ une matrice $m \times n$ et $\mathbf{b}$ un vecteur de $\mathbb{R}^m$ :

$$E[\mathbf{Y}] = A \, E[\mathbf{X}] + \mathbf{b}, \qquad \Sigma_{\mathbf{Y}} = A \, \Sigma_{\mathbf{X}} \, A^T$$

Exemple

Deux capteurs mesurent la même température avec des erreurs $X_1$ et $X_2$, centrées, de variances 9 et 4 (en °C²), et de covariance 2 (ils sont montés côte à côte et subissent en partie les mêmes perturbations). Faut-il faire la moyenne des deux mesures ? L’erreur de la moyenne est $Y = \tfrac{1}{2} X_1 + \tfrac{1}{2} X_2$, soit $A = (\tfrac{1}{2}\ \ \tfrac{1}{2})$ :

$$\mathrm{Var}(Y) = A \, \Sigma \, A^T = \tfrac{1}{4} \times 9 + \tfrac{1}{4} \times 4 + 2 \times \tfrac{1}{4} \times 2 = 2{,}25 + 1 + 1 = 4{,}25$$

La moyenne (4,25) fait à peine mieux que le meilleur capteur seul (4), à cause de la covariance positive. Sans elle, on aurait obtenu 3,25.

Variance de a·X₁ + (1 − a)·X₂ selon le poids a. La covariance positive relève toute la courbe : la moyenne simple (a = 0,5) passe de 3,25 à 4,25, et le meilleur mélange donne plus de poids au capteur le plus précis.

Les vecteurs gaussiens

Un vecteur $\mathbf{X}$ est gaussien si toute combinaison linéaire $a_1 X_1 + \dots + a_n X_n$ de ses composantes suit une loi normale. On note $\mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)$ : sa loi est entièrement déterminée par son vecteur moyen et sa matrice de covariance. Si $\Sigma$ est inversible, sa densité est :

$$f(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} \sqrt{\det \Sigma}} \exp\left(-\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})\right)$$

On reconnaît la densité normale du chapitre 5, où $(x - \mu)^2 / \sigma^2$ est remplacé par la forme quadratique $(\mathbf{x} - \boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})$. Les courbes de niveau de cette densité sont des ellipses (en dimension 2), allongées dans la direction où la variance est la plus grande.

Les ellipses de niveau de N(0, Σ). Leurs axes (flèches) sont les vecteurs propres de Σ, et leurs longueurs sont proportionnelles aux racines des valeurs propres : c'est exactement ce qu'exploite l'ACP.

Les vecteurs gaussiens ont des propriétés remarquables :

• leurs marginales sont gaussiennes ;
• l’image $A\mathbf{X} + \mathbf{b}$ reste gaussienne, de loi $\mathcal{N}(A\boldsymbol{\mu} + \mathbf{b}, A\Sigma A^T)$ ;
• pour un vecteur gaussien, et seulement dans ce cas, composantes non corrélées équivaut à composantes indépendantes : $\Sigma$ diagonale suffit.

Attention

Des marginales toutes gaussiennes ne suffisent pas à faire un vecteur gaussien : il faut que toutes les combinaisons linéaires soient normales. Et la propriété « covariance nulle implique indépendance » ne vaut que pour un vecteur conjointement gaussien.

On simule un vecteur gaussien en R avec la fonction mvrnorm du paquet MASS, et l’on vérifie les formules sur les tirages :

library(MASS)
set.seed(4)
Sigma <- matrix(c(9, 2,
                  2, 4), nrow = 2)            # erreurs des deux capteurs
erreurs <- mvrnorm(100000, mu = c(0, 0), Sigma = Sigma)

cov(erreurs)                                 # proche de Sigma
A <- matrix(c(0.5, 0.5), nrow = 1)
A %*% Sigma %*% t(A)                         # variance théorique de la moyenne : 4.25
var(erreurs %*% t(A))                        # variance observée, proche de 4.25
cor(erreurs)[1, 2]                           # corrélation : 2 / (3 * 2) = 0.33

En Python, numpy.random.default_rng().multivariate_normal(mean, cov, size) joue le même rôle, et numpy.cov estime la matrice de covariance (attention, il attend les variables en lignes par défaut, d’où l’option rowvar=False).

En entretien

« Une covariance nulle ne prouve pas l’indépendance : $Y = X^2$ avec $X$ symétrique en est le contre-exemple type. L’exception, ce sont les vecteurs gaussiens, pour lesquels non corrélé et indépendant sont équivalents. »

En résumé

• La loi conjointe décrit un vecteur aléatoire ; les marginales s’obtiennent en sommant ou en intégrant sur les autres variables, les conditionnelles en normalisant.
• $X$ et $Y$ sont indépendantes si la loi conjointe est le produit des marginales ; la loi d’une somme de variables indépendantes est une convolution.
• $\mathrm{Cov}(X, Y) = E[XY] - E[X]E[Y]$ ; $\rho = \mathrm{Cov}/(\sigma_X \sigma_Y) \in [-1, 1]$ mesure un lien linéaire.
• $\mathrm{Var}(X + Y) = \mathrm{Var}(X) + \mathrm{Var}(Y) + 2\,\mathrm{Cov}(X, Y)$ ; indépendance implique covariance nulle, pas l’inverse.
• La matrice de covariance $\Sigma$ est symétrique semi-définie positive, et $\Sigma_{A\mathbf{X} + \mathbf{b}} = A \Sigma A^T$.
• Un vecteur gaussien $\mathcal{N}(\boldsymbol{\mu}, \Sigma)$ reste gaussien par transformation affine, et ses composantes sont indépendantes si et seulement si $\Sigma$ est diagonale.

Exercices

Exercice · Lire un tableau conjoint

facile

À partir du tableau de la station météo, calculez :

1. $P(Y = 1 \mid X = 0)$ et $P(Y = 1 \mid X = 2)$ ;
2. la probabilité qu’il pleuve au moins un jour et qu’il n’y ait pas de vent fort.

Voir la solution

1. $P(Y = 1 \mid X = 0) = 0{,}05 / 0{,}35 \approx 0{,}14$ et $P(Y = 1 \mid X = 2) = 0{,}20 / 0{,}35 \approx 0{,}57$. Le vent fort est quatre fois plus fréquent les week-ends entièrement pluvieux que les week-ends secs, ce qui confirme la dépendance.
2. On additionne les cases $(X = 1, Y = 0)$ et $(X = 2, Y = 0)$ : $0{,}20 + 0{,}15 = 0{,}35$.

En code

Python

import pandas as pd

# Loi conjointe : lignes X (jours de pluie), colonnes Y (vent fort)
conjointe = pd.DataFrame([[0.30, 0.05],
                          [0.20, 0.10],
                          [0.15, 0.20]], index=[0, 1, 2], columns=[0, 1])

# Loi conditionnelle de Y sachant X : on normalise chaque ligne
cond_y_sachant_x = conjointe.div(conjointe.sum(axis=1), axis=0)
print(cond_y_sachant_x.loc[[0, 2], 1].round(3).tolist())   # [0.143, 0.571]

# Au moins un jour de pluie et pas de vent fort
print(round(conjointe.loc[[1, 2], 0].sum(), 2))           # 0.35

R

# Loi conjointe : lignes X (jours de pluie), colonnes Y (vent fort)
conjointe <- matrix(c(0.30, 0.05,
                      0.20, 0.10,
                      0.15, 0.20), nrow = 3, byrow = TRUE,
                    dimnames = list(X = 0:2, Y = 0:1))

# Loi conditionnelle de Y sachant X : on normalise chaque ligne
cond_y_sachant_x <- prop.table(conjointe, margin = 1)
cond_y_sachant_x[c("0", "2"), "1"]      # 0.143 0.571

# Au moins un jour de pluie et pas de vent fort
sum(conjointe[c("1", "2"), "0"])        # 0.35

Exercice · Variance d'une combinaison

moyen

Un vecteur $(X_1, X_2)$ a pour matrice de covariance $\Sigma = \begin{pmatrix} 2 & -1 \\ -1 & 3 \end{pmatrix}$.

1. Calculez $\mathrm{Var}(3X_1 + 4X_2)$.
2. Calculez $\mathrm{Cov}(X_1 + X_2, X_1 - X_2)$.
3. Calculez le coefficient de corrélation entre $X_1$ et $X_2$.

Voir la solution

1. Avec $A = (3\ \ 4)$ : $\mathrm{Var} = 9 \times 2 + 16 \times 3 + 2 \times 3 \times 4 \times (-1) = 18 + 48 - 24 = 42$.
2. La covariance est bilinéaire : $\mathrm{Cov}(X_1 + X_2, X_1 - X_2) = \mathrm{Var}(X_1) - \mathrm{Cov}(X_1, X_2) + \mathrm{Cov}(X_2, X_1) - \mathrm{Var}(X_2) = 2 - 3 = -1$.
3. $\rho = -1 / \sqrt{2 \times 3} \approx -0{,}41$.

En code

Python

import numpy as np

sigma = np.array([[2, -1],
                  [-1, 3]])

a = np.array([3, 4])
print(a @ sigma @ a)                    # 42   Var(3 X1 + 4 X2)

# Cov(u.X, v.X) = u' Sigma v, avec u = (1, 1) et v = (1, -1)
u, v = np.array([1, 1]), np.array([1, -1])
print(u @ sigma @ v)                    # -1

print(sigma[0, 1] / np.sqrt(sigma[0, 0] * sigma[1, 1]))   # -0.408

R

sigma <- matrix(c(2, -1,
                  -1, 3), nrow = 2)

a <- c(3, 4)
t(a) %*% sigma %*% a                    # 42   Var(3 X1 + 4 X2)

# Cov(u.X, v.X) = u' Sigma v, avec u = (1, 1) et v = (1, -1)
u <- c(1, 1); v <- c(1, -1)
t(u) %*% sigma %*% v                    # -1

cov2cor(sigma)[1, 2]                    # -0.408

Exercice · Somme de deux attentes exponentielles

difficile

À un arrêt, les bus arrivent selon des temps exponentiels de paramètre $\lambda$, indépendants. Vous attendez le deuxième bus (le premier est bondé). Votre attente est $Z = T_1 + T_2$, avec $T_1$ et $T_2$ indépendantes de loi $\mathcal{E}(\lambda)$.

1. Calculez la densité de $Z$ par convolution.
2. Déduisez-en $E[Z]$ et $\mathrm{Var}(Z)$ sans calculer d’intégrale supplémentaire.

Voir la solution

1. Les densités sont nulles sur les négatifs, donc pour $z \geq 0$, l’intégrale se restreint à $0 \leq y \leq z$ :

$$f_Z(z) = \int_0^z \lambda e^{-\lambda(z - y)} \, \lambda e^{-\lambda y} \, dy = \lambda^2 e^{-\lambda z} \int_0^z dy = \lambda^2 \, z \, e^{-\lambda z}$$

C’est la loi gamma de paramètres $(2, \lambda)$. Contrairement à l’exponentielle, sa densité est nulle en 0 : attendre deux bus prend forcément un peu de temps.

2. Par linéarité, $E[Z] = 2/\lambda$. Par indépendance, les variances s’ajoutent : $\mathrm{Var}(Z) = 2/\lambda^2$.

Vérification en R, avec $\lambda = 0{,}1$ par minute : z <- rexp(1e5, 0.1) + rexp(1e5, 0.1), puis mean(z) (proche de 20) et var(z) (proche de 200) ; hist(z, freq = FALSE); curve(dgamma(x, shape = 2, rate = 0.1), add = TRUE) superpose la densité trouvée.

En code

Python

import numpy as np
from scipy import stats

lam = 0.1                                     # un bus toutes les 10 minutes
rng = np.random.default_rng(0)
z = rng.exponential(1 / lam, 100_000) + rng.exponential(1 / lam, 100_000)
print(z.mean(), z.var())                      # 19.97  200.5, proches de 20 et 200

# Comparaison avec la loi gamma(2, lambda) : même proportion sous 15 minutes ?
gamma_2 = stats.gamma(a=2, scale=1 / lam)
print(np.mean(z <= 15), gamma_2.cdf(15))      # 0.4443  0.4422

R

lambda <- 0.1                                 # un bus toutes les 10 minutes
set.seed(0)
z <- rexp(1e5, lambda) + rexp(1e5, lambda)
c(mean(z), var(z))                            # 19.99  200.4, proches de 20 et 200

# Comparaison avec la loi gamma(2, lambda) : même proportion sous 15 minutes ?
c(mean(z <= 15), pgamma(15, shape = 2, rate = lambda))   # 0.4433  0.4422

# Histogramme et densité trouvée par convolution
hist(z, breaks = 60, freq = FALSE)
curve(dgamma(x, shape = 2, rate = lambda), add = TRUE, lwd = 2)