Loi normale: cas multivarié

On considère ici \mathbb{R}^d muni du produit scalaire euclidien \langle \cdot, \cdot \rangle et de la norme euclidienne \|\cdot\| associée.

Rappel sur les vecteurs aléatoires

Si \mathbf{X} = (X_1, \dots, X_d) \in \mathbb{R}^d est un vecteur aléatoire, son espérance est définie coordonnée par coordonnée : \mathbb{E}[\mathbf{X}] = (\mathbb{E}[X_1], \dots, \mathbb{E}[X_d]) \in \mathbb{R}^d\ quantité qui existe dès que chaque espérance est bien définie. De plus, si \mathbb{E}[X_j^2] < \infty pour tout j\in \llbracket 1, d\rrbracket, on peut alors définir les covariances pour tout (i,j) \in \llbracket 1, d\rrbracket^2 : \textrm{cov}(X_i, X_j) = \mathbb{E}[(X_i- \mathbb{E}[X_i]) (X_j - \mathbb{E}[X_j])] \enspace, quantités que l’on rassemble dans une matrice appelée matrice de variance-covariance : \Sigma = (\textrm{cov}(X_i, X_j))_{1 \leq i,j \leq d} \in \mathbb{R}^{d \times d} \,. On peut montrer que cette matrice est symétrique et semi-définie positive. En particulier, si les X_j sont indépendants, alors \Sigma est une matrice diagonale.

Note

Un point important. Si (X,Y) est un vecteur aléatoire, il ne suffit pas de connaître les marginales X et Y pour caractériser entièrement le vecteur. Par exemple, si X et Y suivent toutes les deux une loi normale, alors on peut avoir par exemple X=Y, ou bien X indépendant de Y, et ces deux cas modélisent clairement deux vecteurs aléatoires de lois distinctes.

Vecteurs gaussiens

Définition 1 (Vecteur gaussien) Un vecteur aléatoire \mathbf{X} = (X_1, \dots, X_d)^\top \in \mathbb{R}^d est un vecteur gaussien si pour tout {\alpha} = (\alpha_1, \dots, \alpha_d)^\top, la variable aléatoire réelle \langle {\alpha}, \mathbf{X} \rangle = \alpha_1 X_1 + \cdots + \alpha_d X_d \enspace, suit une loi normale.

En particulier chaque (loi marginale) X_j suit alors une loi gaussienne (choisir ci-dessus \alpha = e_j, les autres égaux à 0).

Cependant, il ne suffit pas que les X_j soient des gaussiennes pour que le vecteur X soit un vecteur gaussien. Par exemple, si X suit une loi normale centrée réduite et \varepsilon une loi uniforme (discrète) sur \{-1,1\}, alors on peut montrer que \varepsilon X suit encore une loi normale centrée réduite. En effet, pour tout t \in \mathbb{R}, on a

\begin{align*} \mathbb{P}(\varepsilon X \leq t) & = \mathbb{P}(X \leq t) \mathbb{P}(\varepsilon = 1) + \mathbb{P}(-X \leq t) \mathbb{P}(\varepsilon = -1)\\ & = \tfrac{1}{2} \mathbb{P}(X \leq t) + \tfrac{1}{2} \mathbb{P}(-X \leq t) = \mathbb{P}(X \leq t) \enspace. \end{align*}

Cependant, X + \varepsilon X prend la valeur 0 avec probabilité 1/2 donc ne suit pas une loi normale. Ainsi le vecteur (X, \varepsilon X)^{\top} n’est pas un vecteur gaussien bi-dimensionnel.

Soit \mathbf{X} un vecteur gaussien. Notons {\mu} son espérance et \Sigma sa matrice de variance-covariance. En reprenant les notations de la définition, la variable aléatoire \langle {\alpha}, \mathbf{X} \rangle vérifie

\begin{align*} \mathbb{E}[\langle {\alpha}, \mathbf{X} \rangle] & = \mathbb{E}[\alpha_1 X_1 + \cdots + \alpha_d X_d] \\ & = \alpha_1 \mathbb{E}[X_1] + \cdots + \alpha_d \mathbb{E}[X_d] & = \langle {\alpha}, {\mu} \rangle\,, \end{align*}

et

\begin{align*} \mathrm{var}(\langle {\alpha}, \mathbf{X} \rangle) & = \mathrm{var}(\alpha_1 X_1 + \cdots + \alpha_d X_d)\\ & = \mathrm{cov}(\alpha_1 X_1 + \cdots + \alpha_d X_d, \alpha_1 X_1 + \cdots + \alpha_d X_d) \\ & = \sum_{1 \leq i,j \leq d} \alpha_i \mathrm{cov}(X_i, X_j) \alpha_j = {\alpha}^\top \Sigma {\alpha} \end{align*}

Ainsi, \langle {\alpha}, \mathbf{X} \rangle \sim \mathcal{N}\left(\langle {\alpha}, {\mu} \rangle, {\alpha}^\top \Sigma {\alpha}\right).

Fonction caractéristique d’un vecteur gaussien

La fonction caractéristique d’un vecteur gaussien d’espérance {\mu} et de matrice de variance-covariance \Sigma est donnée par \phi_\mathbf{X}({\alpha}) = \mathbb{E}[e^{i \langle {\alpha}, \mathbf{X} \rangle}] = \exp\Big(i \langle {\alpha}, {\mu} \rangle - \frac{{\alpha}^\top \Sigma {\alpha}}{2}\Big)\,, \quad {\alpha} \in \mathbb{R}^d\,. où on a utilisé l’expression de la fonction caractéristique d’une variable aléatoire de loi normale \mathcal{N}(\langle {\alpha}, {\mu} \rangle, {\alpha}^\top \Sigma {\alpha}). Ainsi, \phi_\mathbf{X} est entièrement déterminée par les quantités {\mu} et \Sigma. Comme cette fonction caractérise la loi de \mathbf{X}, on en déduit que la loi d’un vecteur gaussien est entièrement caractérisée par son espérance et sa matrice de variance-covariance. On note alors \mathbf{X} \sim \mathcal{N}({\mu}, \Sigma)\,. En particulier, si les variables aléatoires X_1, \dots, X_d sont indépendantes de loi \mathcal{N}(0,1), alors {\mu} = (0,\ldots,0)^\top et \Sigma = \mathrm{Id}_d.

Densité de probabilité

Commençons par le cas {\mu}=0 et \Sigma = \mathrm{Id}_d correspond à la loi gaussienne centrée réduite \mathcal{N}(0, \mathrm{Id}_d). La loi de (X_1,\dots,X_n)^\top correspond alors à la loi produit de n lois gaussiennes centrées réduites indépendantes (pour les gaussiennes la décorrélation implique l’indépendance). Sa densité est donc donnée par \varphi_{0,\mathrm{Id}_d}(x) = \frac{1}{ \sqrt{(2\pi)^d}} \exp\left( -\tfrac{1}{2}x^\top x \right) \enspace.

Proposition 1 (Densité de la loi gaussienne multivariée) Soient {\mu} \in \mathbb{R}^d et \Sigma \in \mathbb{R}^{d \times d} (symétrique et définie positive) et supposons que X \sim \mathcal{N}({\mu},\Sigma). Alors la densité de probabilité de X est donnée pour tout x \in \mathbb{R}^d par

\varphi_{{\mu},\Sigma}(x) = \frac{1}{ \sqrt{(2\pi)^d |\det(\Sigma)|}} \exp\Big( -\tfrac{1}{2}(x-{\mu})^\top\Sigma^{-1}(x - {\mu}) \Big) \enspace.

Preuve. Prenons une matrice L\in \mathbb{R}^{d \times d} telle que LL^\top = \Sigma (par exemple, la décomposition de Cholesky, que nous reverrons ci-dessous). Calculons alors la loi de \mathbf{Y} = \psi(X) \triangleq L \mathbf{X} + {\mu}, pour X\sim \mathcal{N}(0,\mathrm{Id_d}). Pour appliquer la formule du changement de variable nous donne, on calcule \psi^{-1} et son jacobien: \psi^{-1}(y) = L^{-1}(y-{\mu}), et \det(J_{\psi^{-1}}) = \det(L^{-1}) = \det(L)^{-1} = |\det(\Sigma)|^{-1/2}.

En notant que LL^\top \left(L^{-1}\right)^\top L^{-1}=\mathrm{Id}_d, et donc que \left(L^{-1}\right)^\top L^{-1}=\Sigma^{-1}, on en déduit que la densité de \mathbf{Y} est donnée par

\begin{align*} \varphi_{{\mu},\Sigma}(y) & = \varphi_{0,\mathrm{Id}_d}(\psi^{-1}(y)) |\det(J_{\psi^{-1}})| \\ & = \frac{|\det(\Sigma)|^{-1/2}}{ \sqrt{(2\pi)^d }} \exp\left( -\tfrac{1}{2}(y-{\mu})^\top \left(L^{-1}\right)^\top L^{-1}(y - {\mu})\right)\\ & = \frac{1}{ \sqrt{(2\pi)^d |\det(\Sigma)|}} \exp\Big( -\tfrac{1}{2}(y-{\mu})^\top\Sigma^{-1}(y - {\mu}) \Big) \enspace. \end{align*}

Proposition 2 (Transformation affine de vecteurs gaussiens) Soit \mathbf{X} \sim \mathcal{N}({\mu}, \Sigma) un vecteur gaussien sur \mathbb{R}^d, \Omega \in \mathbb{R}^{d' \times d} et {\nu}\in \mathbb{R}^{d'}. Alors, le vecteur aléatoire \mathbf{Y} = \Omega \mathbf{X} + {\nu} est un vecteur gaussien vérifiant \mathbf{Y} \sim \mathcal{N}(\Omega {{\mu}} + {\nu}, \Omega \Sigma \Omega^\top)\,.

Cette proposition se prouve sans peine en utilisant la fonction caractéristique On retrouve en particulier la stabilité par transformation affine établie en dimension 1.

La factorisation de Cholesky

Pour rappel, la factorisation de Cholesky¹ d’une matrice symétrique définie positive est donnée ci-dessous.

¹ André-Louis Cholesky, dit René: (1875-1918) ingénieur topographe et géodésien dans l’armée française, mort des suites de blessures reçues au champs de bataille.

Théorème 1 (Factorisation de Cholesky) Soit \Sigma \in \mathbb{R}^{d \times d} une matrice symétrique définie positive. Alors il existe une matrice triangulaire inférieure L \in \mathbb{R}^{d \times d} telle que \Sigma = LL^\top. La décomposition est unique si l’on impose que les éléments diagonaux de L soient strictement positifs.

Preuve. La factorisation de Cholesky est une conséquence directe de la méthode du pivot de Gauss. Le détail est donné par exemple dans (Th. 4.4.1, Ciarlet 2006).

Ciarlet, P. G. 2006. Introduction à l’analyse numérique matricielle et à l’optimisation. Cours et exercices corrigés. Dunod.

Simulation de vecteurs gaussiens

La simulation d’un vecteur gaussien de paramètres {\mu} = (0,\ldots,0)^\top et \Sigma = \mathrm{Id}_d ne pose pas de problème : il suffit de simuler X_1,\dots, X_d, d variables aléatoires indépendantes de loi normale centrée réduite. En effet, le vecteur \mathbf{X} = (X_1,\dots, X_d)^\top est alors un vecteur gaussien de loi \mathcal{N}(0, \mathrm{Id}_d).

Supposons maintenant que l’on veuille simuler un vecteur gaussien de loi \mathcal{N}({\mu}, \Sigma) dans \mathbb{R}^d, {\mu} et \Sigma symétrique définie positive donnés.

Approche par la factorisation de Cholesky

La matrice \Sigma étant symétrique, elle peut s’écrire comme \Sigma = LL^\top où L est une matrice triangulaire inférieure de taille d \times d. Grâce à la décomposition de Cholevsky et en reprenant les éléments de la preuve de la Proposition 1, on peut écrire \mathbf{Y} = L \mathbf{X} + {\mu} où \mathbf{X} \sim \mathcal{N}(0, \mathrm{Id}_d) et vérifier que \mathbf{Y} \sim \mathcal{N}({\mu}, \Sigma).

Approche par la décomposition spectrale de \Sigma

La matrice \Sigma étant symétrique, elle se diagonalise en base orthonormée : il existe une matrice orthogonale P telle que \Sigma = P \mathrm{diag}(\lambda_1 \ldots, \lambda_d) P^{-1} = P \mathrm{diag}(\lambda_1 \ldots, \lambda_d) P^\top\,, où \lambda_1, \ldots, \lambda_d \geq 0 sont les valeurs propres de \Sigma qui est semi-définie positive. On pose alors R = P \mathrm{diag}(\sqrt \lambda_1 \ldots, \sqrt \lambda_d) qui est une racine carrée matricielle de \Sigma au sens où \Sigma = R R ^\top. On part alors d’un vecteur gaussien centrée réduit \mathbf{X}_0 \sim \mathcal{N}(0, \mathrm{Id}_d) que l’on sait simuler (par exemple avec la méthode de Box-Müller). La proposition Proposition 2 assure alors que le vecteur \mathbf{X} = R \mathbf{X}_0 + {\mu} est un vecteur gaussien de loi \mathcal{N}({\mu}, \Sigma).

Vecteurs gaussiens : cas bidimensionnel

En dimension p=2, la matrice de covariance \Sigma peut toujours s’écrire comme suit, et la visualisation suivante montre l’impact des différents paramètres sur la densité de probabilité. \Sigma = \begin{pmatrix}\cos(\theta) & - \sin(\theta)\\ \sin(\theta)& \cos(\theta) \end{pmatrix} \cdot \begin{pmatrix}\sigma_1 & 0\\ 0 & \sigma_2 \end{pmatrix}\cdot \begin{pmatrix} \cos(\theta) &\sin(\theta)\\ -\sin(\theta)& \cos(\theta)\end{pmatrix} où \theta est l’angle de rotation et \sigma_1 et \sigma_2 les écarts-types des marginales (dans le repère orthonormal après rotation).

Plotly = require('plotly.js-dist');
dists = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists@umd/browser.js' );
jstat = require('jstat');
math = require("mathjs");
// see source here: https://github.com/stdlib-js/stats-base-dists/tree/umd
// continuous case

viewof inputs = Inputs.form([
  Inputs.range([-1, 1], {label: tex`\mu_1`, step: 0.1}, {value: 0}),
  Inputs.range([-1, 1], {label: tex`\mu_2`, step: 0.1}, {value: 0}),
  Inputs.range([0.1, 2], {label: tex`\sigma_1`, step: 0.1, value: 1}),
  Inputs.range([0.1, 2], {label: tex`\sigma_2`, step: 0.1, value: 1}),
  Inputs.range([0, 6.29], {label: tex`\theta`, step: 0.01, value: 0}),
  Inputs.range([0, 1000], {label: tex`n`, step: 1, value: 10}),
  Inputs.button("Re-Tirage")
])


mu1 = inputs[0];
mu2 = inputs[1];
sigma1 = inputs[2];
sigma2 = inputs[3];
theta = inputs[4];
n_samples = inputs[5];


function create_sigma(theta, sigma1, sigma2){
  const mat_rot = math.matrix([[math.cos(theta), -math.sin(theta)], [math.sin(theta), math.cos(theta)]]);
  const mat_sigma = math.matrix([[sigma1**2, 0], [0, sigma2**2]]);
  return math.multiply(mat_rot, math.multiply(mat_sigma, math.transpose(mat_rot)));
}


function mvnpdf(x, mu, Sigma){
  const p = 2;
    return (2*math.pi)**(-p/2)*math.det(Sigma)**(-0.5)*
      math.exp(-0.5*math.multiply(math.multiply( math.transpose(math.subtract(x,mu)), math.inv(Sigma)), math.subtract(x, mu)));
}

{


function normal_rng(mu, Sigma, n=100) {
    // Compute the Cholesky decomposition of Sigma
    var cholesky = jstat.cholesky(Sigma);

    // Generate the samples
    var samples = Array.from({length: n}, () => {
        var x = jstat.normal.sample(0, 1);
        var y = jstat.normal.sample(0, 1);
        // Transform the standard normal random variables using the Cholesky decomposition
        var transformedX = mu[0] + cholesky[0][0] * x + cholesky[0][1] * y;
        var transformedY = mu[1] + cholesky[1][0] * x + cholesky[1][1] * y;
        return [transformedX, transformedY];
    });   

    return samples;
}
var Sigma = create_sigma(theta, sigma1, sigma2).toArray();
var mu = [mu1, mu2];

var samples = normal_rng(mu, Sigma, 1000);
var npoints=100, mini = -5, maxi=5, x = new Array(npoints), y = new Array(npoints), z = new Array(npoints), i, j;

//  Densité:
for(var i = 0; i < npoints; i++) {
    x[i] = mini + i * (maxi - mini) / (npoints - 1);
    y[i] = mini + i * (maxi - mini) / (npoints - 1);
    z[i] = new Array(npoints);
    }

for(var i = 0; i < npoints; i++) {
    for(j = 0; j < npoints; j++) {

        z[j][i] = mvnpdf([x[i], x[j]], mu, Sigma);
    }

}
  

{

var trace1 = {
        x: samples.slice(0, n_samples).map(sample => sample[0]),
        y: samples.slice(0, n_samples).map(sample => sample[1]),
        mode: 'markers',
        type: 'scatter',
        marker: {
            color: 'rgba(0,0,0,0.5)',
            size: 5,
        },

        xaxis: 'x2',
}

var trace22 = {
        x: x,
        y: y,
        z: z,
        type: 'surface',
        colorscale: 'Oranges',
        showscale: false,
        color: {
            legend: false,
            label: "pdf",
        },

}

var trace21 = {
        x: x,
        y: y,
        z: z,
        type: 'contour',
        colorscale: 'Oranges',
        color: {
            legend: true,
            label: "pdf",
        },
        blur: 4,
        xlim: [-5, 5],
        ylim: [-5, 5],
        xaxis: 'x2',
}


var data = [
  trace21,
  trace22,
  trace1,
  ];


  var layout = {
       yaxis2: {
          range: [-5, 5],
          autorange: false,
        },
       xaxis2: {
          range: [-5, 5],
          autorange: false,
        },

      scene: {
          camera: {
              eye: {
                  x: 0.5,
                  y: 1.2,
                  z: 1.5,
              }
          }
      },
    grid: {
      rows: 1,
      columns: 2,
      subplots: [['xy','x2y']],
    },
    showlegend: false,

  };

    var config = {responsive: true}
    const div = DOM.element('div');
    Plotly.newPlot(div, data, layout, config);
    return div;
  }

}