Génération de variables aléatoires à distribution gaussienne

Lorsque l’on développe une application Web ou mobile, il n’est pas rare d’avoir besoin de générer des valeurs aléatoires. C’est un fait. Mais dans nombre de situations, la distribution uniforme offerte par les fonctions natives – et c’est le cas de Math.random() en Javascript – ne répond pas aux besoins. En d’autres termes, tous les tirages sur l’intervalle cible ( [0,1[ pour Math.random() par exemple) sont équiprobables, et cela ne permet parfois pas de générer un ensemble de données plausibles. C’est par exemple le cas lorsque l’on veut simuler des phénomènes naturels. Voyons donc comment implémenter une fonction de génération de nombres aléatoires suivant une loi normale, ou Gaussienne. En Javascript. Parce que pourquoi? Parce que parce que.

Distribution uniforme

En générant un nombre aléatoire avec la fonction Math.random(), j’ai donc autant de chance de tirer n’importe quelle valeur sur l’intervalle [0,1[. « Prouve-le! ». « D’accord ».

Commençons par ré-expliquer un peu ce qu’est une distribution uniforme, car on me souffle dans l’oreillette que cela n’est peut-être pas clair pour tout le monde. Pour faire simple, cela signifie que les probabilités de tirages de toutes les valeurs de l’intervalle cible sont égales. Prenons l’exemple d’un dé pour illustrer cela. Sur un tirage, j’ai autant de chances d’obtenir un 1, un 2, un 3, un 4, un 5, ou un 6. Soit P(1) = P(2) = P(3) = P(4) = P(5) = P(6) = 1/6. De fait, si je procède à un nombre suffisamment grand de tirages avec ce même dé, j’aurai, à terme, obtenu approximativement obtenu autant de 1, que de 2, que de 3… la suite vous la connaissez. Si alors j’illustre par un histogramme le nombre de tirages pour chaque face, tous mes petits bâtons auront la même taille, à peu de chose près, et c’est cela une distribution uniforme.

Écrivons donc un petit algo – eeeextrêmement compliqué, ça frôle le génie. je sais. – générant un tableau de, disons, 100000 valeurs aléatoires avec cette fonction.

var array = [],
    length = 1e5;

function createUniformRandomArray() {
    // create an array of 'length' elements
    for ( var i = 0; i < length; i++ ) {
        // filled with random values generated by Math.random() function
        array.push( Math.random() );
    }
}

createUniformRandomArray();

Cela fait, visualisons la distribution, i.e. la densité des tirages. Pour ce faire, nous allons utiliser la librairie Plotly (https://plot.ly/). L’utilisation de celle-ci n’étant pas le sujet de ce billet je ne rentrerai pas dans le détail, mais, comme vous pourrez le constater, c’est assez trivial.

// create density data for graph
var density_data = {
    x: array,
    name: 'x density',
    type: 'histogram'
};

// plot the density distribution using Plotly
Plotly.newPlot('uniform-distribution-graph', [ density_data ] );

Et voici donc la fameuse représentation graphique de la distribution pour l’un de mes tirages:

Plutôt uniforme n’est-ce pas? Eh bien ça alors, on ne s’y attendait pas….

Si vous ne me croyez pas sur parole – mécréants – vous pouvez jouer un peu sur le CodePen suivant:

SimplX – Uniform random distribution

Nous devrions donc être désormais relativement d’accord sur le fait que Math.random() génère un ensemble de données aléatoires avec une répartition uniforme. Voyons maintenant comment nous pouvons utiliser cette même fonction pour tirer des valeurs suivant une loi normale, ou Gaussienne, afin de « centrer » la répartition sur une valeur moyenne.

Distribution  Gaussienne

Pour arriver à nos fins, nous allons nous reposer sur la transformation de Box-Muller, mais notons qu’il ne s’agit pas de la seule méthode, bien que très certainement la plus utilisée.

Cette méthode permet, à partir de deux nombres aléatoires indépendants, prenant valeur dans [0,1] avec une distribution uniforme, de générer deux variables aléatoires, indépendantes elles aussi, suivant une distribution Gaussienne centrée sur 0, avec une déviation standard de 1. Ici encore, je ne rentrerai pas dans le détail mathématique de tout cela, ça n’est pas le sujet, mais soyez rassurés si vous êtes un peu perdus, la suite devrait très largement vous éclairer.

Dans sa forme générale, voici à quoi ressemble cette petite bête:

où X1 et X2 sont les deux variables indépendantes générées, suivant une distribution gaussienne, et r1 et r2 les deux valeurs aléatoires initiales, avec distribution uniforme donc. Suivez un peu.

NdlR: Il est bien sûr possible de n’utiliser qu’une seule de ces deux variables, et c’est d’ailleurs ce que nous allons faire ici, mais il est toujours bon de savoir qu’on peut en avoir une seconde pour le même prix, et qu’elles sont, je le répète, totalement indépendantes.

Pour des raisons de performance et de stabilité dans certaines situations, ce n’est pas cette forme que nous allons utiliser, mais sa forme polaire, dont l’implémentation est la suivante:

function gaussianRandom() {
    var r1, r2, w, X1, X2;

    do {
        r1 = 2 * Math.random() - 1;
        r2 = 2 * Math.random() - 1;
        w = r1 * r1 + r2 * r2;
    } while ( w >= 1 );

    w = Math.sqrt( ( -2 * Math.log( w ) ) / w );

    X1 = r1 * w;
    X2 = r2 * w;

    return X1;
}

Pour ceux qui souhaiteraient se rafraîchir la mémoire sur la forme polaire des nombres complexes et la manière de passer de la forme cartésienne à la forme polaire, voici un petit article Wikipedia. Cela n’est pas indispensable pour comprendre la suite, mais fera sûrement le bonheur des plus curieux!

Maintenant, nous pouvons faire la même chose que dans la précédente section, générer un tableau de ces valeurs, puis afficher leur distribution.

var array2 = [];

function createGaussianRandomArray() {
    // reset the array
    array2 = [];

    // fill the array with 'length' elements generated by
    // the previously created gaussianRandom() function
    for ( var i = 0; i < length; i++ ) {
        array2.push( gaussianRandom() );
    }
}

createGaussianRandomArray();

// create density data for graph
var density_data = {
    x: array2,
    name: 'x density',
    type: 'histogram'
};

// plot the density distribution using Plotly
Plotly.newPlot('gaussian-distribution-graph', [ density_data ] );

Et nous constatons avec un immense plaisir que l’on obtient une jolie distribution Gaussienne. C’est beau.

Évidemment, les plus pointilleux me rétorquerons à l’unisson « c’est bien beau tout ça, mais visiblement je tire des valeurs entre « -4 et des poussières », et « +4 et des poussières », et non sur [0,1]…! ». Certes. Alors pour les plus exigeants, appliquons une dernière transformation. En même temps je vous avais promis une moyenne à 0 et une déviation standard de 1, mais bon… ok. Une simple normalisation devrait très certainement faire l’affaire, non?

Mais avant de s’attaquer à ce sujet, pour les plus sceptiques, une petite explication sur ces fameux tirages sur [-4, +4]. Nous avons donc un algorithme qui génère des tirages suivant une distribution gaussienne avec une déviation standard de 1. Alors pourquoi ce « -4 et des poussières », et « +4 et des poussières »…? Pour le comprendre, rappelons les percentiles d’une distribution normale. Un joli dessin valant toujours mieux que pléthore de vilains mots, voyons cela en image:

Notre déviation standard, ou sigma, étant de 1, nous aurons donc 99,7% de tirages entre -3 et +3, mais quelques vilains petits canards (0.3%) au delà… Donc entre -4 et +4, voir encore au delà, on n’est jamais à l’abri…

Un peu plus clair? Passons donc à cette sombre histoire de normalisation…

Pour faire simple, l’idée est de passer de valeurs comprises dans un intervalle [a,b] à des valeurs comprises dans un intervalle [0,1]. Nous commençons donc par décaler l’ensemble des valeurs dans l’intervalle [0, b-a]. Cela revient simplement à effectuer une translation de nos valeurs, soit à soustraire a à toutes les valeurs. Puis il ne reste plus qu’à « compresser » tout cela pour ramener l’ensemble des valeurs dans [0,1], et il suffit désormais pour cela simplement de diviser toutes les valeurs par la distance entre a et b, soit (b-a).

Rajoutons que pour être totalement rigoureux, nous allons également modifier les bornes de l’intervalle de départ pour qu’il soit bien centré en 0. Pour ce faire, nous déterminons d’abord en valeur absolue la plus grande valeur, de min ou de max, que nous appellerons boundary, et appliquons la normalisation sur [-boundary, +boundary].

function normalize( array ) {
 
 var min = Math.min.apply( null, array ),
 max = Math.max.apply( null, array );
 
 var boundary = Math.max(-min, max);
 
 min = -boundary;
 max = boundary;

 var new_array = [];

 for ( var i = 0; i < array.length; i++ ) {
   new_array.push( ( array[i] - min ) / ( max - min ) );
 }

 return new_array;
}

Puis on représente la distribution du fraîchement créé tableau normalisé.

var normalized_array = normalize( array2 );
    
// create density data for graph
var density_data = {
    x: normalized_array,
    name: 'x density',
    type: 'histogram'
};
    
// plot the density distribution using Plotly
Plotly.newPlot('gaussian-distribution-graph', [ density_data ] );

Comme diraient nos amis américains: « et voilà! ». Et nous aussi, d’ailleurs. Voleurs.

Une belle distribution en loi normale encore, mais prenant valeur dans [0,1], et centrée sur 0.5. Et là on est bien.

Ici encore, je vous invite à constater tout cela de vos propres yeux sur le CodePen, et à vous amuser un peu avec le code. Ça ne mange pas de pain, et ça fait toujours plaisir.

SimplX – Gaussian random distribution

BONUS: le théorème central limite

Le théorème central limite stipule que toute somme de variables aléatoires à distribution uniforme tend vers une variable aléatoire à distribution gaussienne. Sans rentrer dans le détail de ce théorème, nous voyons donc ici une autre manière – bien que plus gourmande en calcul – de générer de telles variables.

Je terminerai donc cet article en vous laissant avec un dernier CodePen illustrant ce théorème, en guise de cerise sur le gâteau… gaussien.