Génération de nombre aléatoire suivant une loi gaussienne

**millie** · 08/05/2007, 10h26

Bonjour,

Il peut parfois arriver que l'on ait besoin de générer aléatoirement des nombres aléatoires suivant une loi gaussienne, je fais donner l'algorithme basé sur la technique de Box-Muller.

Tout d'abord, on suppose que l'on dispose d'une fonction Aleatoire0-1 permettant de générer aléatoirement des réels entre 0 et 1 strictement. (par exemple Math.random de Java si je ne m'abuse).

Si vous ne disposez pas d'une telle fonction, alors vous disposez surement d'une fonction (du genre rand) qui permet de générer un entier aléatoire quelconque.
Pour simuler la fonction Aleatoire0-1, il est possible d'écrire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
fonction Aleatoire0-1() -> réel
{
  /*on choisi un entier entre 1 et un nombre très grand (par exemple 65001)
  Entier i = rand() % 65000+1; 
 /*on divise par 65001 pour avoir un réel entre 0 et 1 strictement*/
  Reel f = i / 65001; 
  retourner f
}

A présent, on peut générer un nombre aléatoire suivant une loi gaussienne d'espérance mean et de variance sigma²

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
fonction gaussienne(Réel sigma, Réel mean)
 
 Réel u1 = Aleatoire0-1();
 Réel u2 = Aleatoire0-1();
 
 retourner(racine(-2 * ln(u1)) * cos( 2 * pi * u2) * sigma + mean
 
}

L'implémentation de Java propose le code suivant (on en parle après) (code sous licence GPL) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
    synchronized public double nextGaussian() {
        // See Knuth, ACP, Section 3.4.1 Algorithm C.
        if (haveNextNextGaussian) {
    	    haveNextNextGaussian = false;
    	    return nextNextGaussian;
    	} else {
            double v1, v2, s;
    	    do {
                v1 = 2 * nextDouble() - 1; // between -1 and 1
            	v2 = 2 * nextDouble() - 1; // between -1 and 1
                s = v1 * v1 + v2 * v2;
    	    } while (s >= 1 || s == 0);
    	    double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
    	    nextNextGaussian = v2 * multiplier;
    	    haveNextNextGaussian = true;
    	    return v1 * multiplier;
        }
    }

**pseudocode** · 08/05/2007, 13h06

Pour les java-fan-boys, la classe Random permet de generer une valeur aléatoire gaussienne de moyenne (mean) 0 et d'ecart-type (sigma) 1.0

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
Random generator;
generator = new Random();
double num = generator.nextGaussian();

Ensuite, une simple transformation affine permet de modifier la moyenne et l'ecart-type:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

double num2 = ( num  * sigma ) + mean

**PRomu@ld** · 09/05/2007, 07h52

J'ai l'habitude d'utiliser Mersenne Twister :

http://www.math.sci.hiroshima-u.ac.j...eversions.html

**Jean-Marc.Bourguet** · 09/05/2007, 09h07

Promu@Id, l'algo de Mersenne-Twister permet de generer une sequence uniformement repartie, le sujet ici est d'en faire une sequence repartie suivant une loi gaussienne.

millie, j'ai une question et une suggestion. La question, est-ce que les bornes 0 et 1 doivent pouvoir etre generee (auquel cas, il vaut mieux preciser la maniere doit ln(0) doit etre gere)?

La suggestion, enlever la constante 65001 et la remplacer par RAND_MAX ou RAND_MAX+1 suivant la reponse a la question precedente. Utiliser autre chose qu'un diviseur de RAND_MAX+1 introduit un biais d'autant plus important que ce nombre est grand (c'est particulierement visible quand le nombre est plus grand, comme ce le sera pour les distraits utilisant une implementation de C ou RAND_MAX vaut 32767)

**millie** · 09/05/2007, 10h07

Envoyé par Jean-Marc.Bourguet

millie, j'ai une question et une suggestion. La question, est-ce que les bornes 0 et 1 doivent pouvoir etre generee (auquel cas, il vaut mieux preciser la maniere doit ln(0) doit etre gere)?

La suggestion, enlever la constante 65001 et la remplacer par RAND_MAX ou RAND_MAX+1 suivant la reponse a la question precedente. Utiliser autre chose qu'un diviseur de RAND_MAX+1 introduit un biais d'autant plus important que ce nombre est grand (c'est particulierement visible quand le nombre est plus grand, comme ce le sera pour les distraits utilisant une implementation de C ou RAND_MAX vaut 32767)

Pour les bornes ]0,1[. J'avais simplement vu ça sur wikipedia... Maintenant, avoir les bornes ]0,1] ne doit pas changer grand chose. Mais en théorie, on est censé travailler avec des nombres réels, donc la probabilité d'obtenir 0 et 1 est censée être nulle. (juste en théorie, donc il faut quand même règler ce problème là)

L'utilisation de RAND_MAX doit effectivement être beaucoup plus efficace. Je ne connaissais pas, je l'utiliserai donc

Je suis également tombé sur d'autres algorithmes à cette page :
http://www.taygeta.com/random/gaussian.html (il génère deux nombres aléatoires d'un coup). Mais il faut gérer un cas particulier qui n'est pas décrit dans l'algorithme.

**pseudocode** · 09/05/2007, 10h50

Envoyé par millie

Je suis également tombé sur d'autres algorithmes à cette page :
http://www.taygeta.com/random/gaussian.html (il génère deux nombres aléatoires d'un coup). Mais il faut gérer un cas particulier qui n'est pas décrit dans l'algorithme.

C'est le meme algo utilisé par java: http://java.sun.com/j2se/1.4.2/docs/...nextGaussian()

**millie** · 09/05/2007, 11h27

Cet algorithme est peut être moins en carton que le mien

Je ferais quelques tests pour voir la repartition...

Par contre, il est possible de sortir de la boucle avec s =0, et ils font : Math.log(s)/s, et ils ne font pas de cas particulir juste pour s = 0. A moins que java sache gérer la forme particulière ln x /x en 0

**pseudocode** · 09/05/2007, 11h34

Envoyé par millie

Par contre, il est possible de sortir de la boucle avec s =0

Non, on ne peut pas sortir avec s =0:

Envoyé par Random.nextGaussian

do { ...} while (s >= 1 || s == 0);

Envoyé par millie

A moins que java sache gérer la forme particulière ln x /x en 0

Heu... non. Java c'est pas matlab

Au mieux la valeur sera Nan (not a number)

**millie** · 09/05/2007, 12h29

Ah oui quel boulet

. Je sais pas pourquoi je pensais le contraire.

**millie** · 09/05/2007, 22h02

Donc, avec ma méthode de base, voici deux exemples de répartitions, pardonnez moi l'ASCIIART

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *
                   *           *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   **    *     *
                   **    *     *
                   **    *     *
                   **    *     *
                   ***   *     *
                   ***   *     *
                   ***   *     *
                   ***   *   * *
                   ***   *  **** *
                   ***   *  ******
                   ***   *  ******
                   ***   *  ******
                   ****  * *******
                   ****  * *******
                   ***** * *******
                  ******** *******
                  ****************
                  ****************
                  ****************
                *******************
                *******************
                *******************
                *******************
                *******************
               *********************
               *********************
               *********************
              **********************
              ***********************
          *  *************************
          ** *************************
       ** ***************************** * *
 
 
 
 
 
 
                         *
                   *     *
                   *     *
                   *     *
                   *     *
                   *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *     *
                   *     *    **
                   *     *  * **
                   **    *  * **
                   ***   *  * **
                   ***   *  * **
                   ***   *  * ** *
                   ***   ** * ** *
                   ****  ** * ** *
                   ****  ** * ** *
                 ****** *** * ** *
                 ****** ***** ** *
                 ****** ***** ** *
                 ****** ******** *
                 ****** ******** *
                 ****** ******** *
                 *****************
                 *****************
                ******************
                *******************
                *******************
                *******************
                ********************
                ********************
                ********************
                ********************  *
               *********************  *
               ********************** *
           *  *********************** *
          *****************************
          ***************************** *

Ce qui n'est pas vraiment terrible.

En utilisant le conseil de Jean Marc (à savoir utiliser RAND_MAX (rand()/(RAND_MAX+1.0f) pour générer un flottant entre 0 et 1)

On obtient ce genre de répartition :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
                        *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                        **
                       *** *
                       *****
                      ******
                      ******
                      ******
                      ******
                      ******
                      *******
                     ******** *
                     ******** *
                     **********
                    ***********
                    ************
                    ************
                    ************
                    ************
                    ************
                    ************
                    ************
                    ************
                   *************
                   *************
                   **************
                  ***************
                  ***************
                  ***************
                  ***************
                  **************** *
                  **************** *
                 ***************** *
                ********************
                ********************
              * ********************
              **********************
          *   ***********************
          **  ***********************

C'est déjà plus sympa.

Et en utilisant la méthode Polaire de Box-Muller implémenté par Java, on obtient ce genre de répartition :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                        ***
                        ***
                        ***
                        ***
                        ***
                        ***
                        ***
                        ***
                        ***
                        *****
                        *****
                        *****
                        *****
                        *****
                    * * *****
                    *********
                    *********
                    *********
                    *********
                    *********
                    *********
                    ********* *
                    ********* *
                   ************
                   ************
                   ************
                   *************
                   *************
                  **************
                  ***************
                  ***************
                  ***************
                  ***************
                ******************
                ******************
                ******************
                ******************
                ******************
                ******************
                ******************
              * *******************
             ** *******************
             ***********************  *
            ***************************
 
 
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         *
                         **
                         **
                         **
                         **
                         **
                         ***
                         ***
                         ***
                         ***
                         ***
                         ***
                         ***
                         ***
                         ***
                        ****
                        ****
                        ****
                       *****
                       *******
                       *******
                     * *******
                     * ******* *
                     * ******* *
                     * ******* *
                     * *********
                     ***********
                   * ***********
                   *************
                   *************
                   *************
                   *************
                 * *************
                 * *************
                 * **************
                 * **************
                 * **************
                 ****************
                 ****************
                 *****************
                 *****************
                 *****************
              *  ******************
              * ********************
              * *********************
            **************************
         *  **************************** *

Donc voilà

**harek** · 19/01/2012, 12h56

j'ai une question:
comment générer des distrubutions non uniformes , sachant que le systéme ne fournit que des générateurs unifrmes ?

Génération de nombre aléatoire suivant une loi gaussienne

Contribuez

Discussions similaires

Partager

Partager