Faire de l'aléatoire convenablement

Dans un intervalle

L'auteur de ce sujet a trouvé une solution à son problème.
Auteur du sujet

Bonjour à tous !

Dans le cadre d'un programme (que je dois absolument coder en c), j'aimerai reproduire le comportement de la fonction random.randrange(), qui permet de retourner un nombre $x$ tel que $ x \in [\text{min}; \text{max}]$. Ça veux dire aussi que $x \in \mathbb{R}$.

Une solution comme suis n'est donc pas envisageable :

1
2
3
double randrange(double min, double max) {
    return min + rand() % (max-min);
}

Déjà parce que le modulo donne un nombre entier, ensuite parce qu'il est dit un peu partout qu'il ne faut pas faire comme ça (j'ai pas trop pigé pourquoi). Du coup, sur les conseils de cette page, j'ai écris ceci :

1
2
3
4
5
6
7
8
9
double randrange(double M, double N) {
    if (M > N) {
        double T = M;
        M = N;
        N = T;
    }

    return M + rand()/(RAND_MAX + 1.) * (N-M);
}

Mais je ne suis pas sur qu'elle soie meilleure (déjà, pourquoi + 1. ?!?). En plus, pour l'application que j'écris, j'ai besoin de "vrai" aléatoire ou de ce qui s'en rapproche de plus (donc pas une distribution gaussienne), et la plupart des messages sur internet parlent peu ou prou de "casser l'aléatoire" quand on écrit une fonction du genre.

J'ai aussi lu à propos de drand48(), mais je suis pas sous OpenBSD.

Est ce que ma fonction est correcte ? Est ce qu'il y a des manières plus intelligentes de faire ? D'avance merci :)

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0
Staff

Comme je sais que la fonction random() de python est une fonction uniforme, j'ai donc été voir comment ils faisaient et voici ce que me dit le code de cpython :

 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
/* generates a random number on [0,0xffffffff]-interval */

static PY_UINT32_T
genrand_int32(RandomObject *self)
{
    PY_UINT32_T y;
    static const PY_UINT32_T mag01[2] = {0x0U, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */
    PY_UINT32_T *mt;
    mt = self->state;
    if (self->index >= N) { /* generate N words at one time */
        int kk;

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];

        self->index = 0;
    }

    y = mt[self->index++];
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680U;
    y ^= (y << 15) & 0xefc60000U;
    y ^= (y >> 18);
    return y;
}

j'espère que ça t'intéresse.

+0 -0
Staff

Cette réponse a aidé l'auteur du sujet

Pour le côté correcte et programmation en C, je ne peux pas dire, pour le côté bon aléatoire, ça dépend du nombre de tirage que tu veux faire. Si c'est une centaine, tu peux partir sur la fonction rand, si c'est plusieurs de millions (algorithme de type Monte-Carlo, par exemple), il ne faut surtout pas. Au bout d'un moment, les nombre généré par rand ne sont plus correctement distribués. Il faut alors se tourner vers d'autres générateurs de nombres aléatoires, comme Mersenne Twister.

Le rand de python ou le module random de C++ utilisent Mersenne. Pas le rand du C.

Édité par Gabbro

Hier, dans le parc, j'ai vu une petite vieille entourée de dinosaures aviens. Je donne pas cher de sa peau.

+1 -0

Si le but est de ne pas se prendre trop le chou et que tu es sous un système Unix, lire le fichier /dev/random devrait suffire.

Amateur de Logiciel Libre et de la distribution GNU/Linux Fedora.

+0 -0
Auteur du sujet

J'ai écris ça un peu vite ce matin, je vois :p

Déjà, ce que je craignais est arrivé :

Si c'est une centaine, tu peux partir sur la fonction rand, si c'est plusieurs de millions (algorithme de type Monte-Carlo, par exemple), il ne faut surtout pas. Au bout d'un moment, les nombre généré par rand ne sont plus correctement distribués.

Gabbro

C'est bel et bien un algo de type Monte Carlo que j'essaye d'écrire. Je me doutais bien qu'il y avait un truc du genre qui allait sortir, mais les papiers qui décrivent Monte Carlo partent du principe que tu a un Bac+18 en statistique et que ce genre de chose est obvious. Merci, donc :)

@Kje et @artragis : merci. Sauf que je me suis rendu compte que randrange() reste dans le domaine des entiers, or je veux un réel ;) Après, le code de génération d'entier de cpython m'as l'air de correspondre à ce que Gabbro mentionnait comme étant "Mersenne Twister", donc je peux éventuellement leur emprunter … Mais faut que je sois sur que ma fonction randrange() du dessus me détruise pas la qualité de l'aléatoire généré.

@Renault: de fait, et à priori, je resterai sur de l'UNIX tout du long. Néanmoins, on ne sais jamais (question de rapidité, /dev/urandom est pas "lent" vu qu'il y a une lecture à faire ?).

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

Je ne comprens pas très bien ton message : a priori, ce que tu demande correspond à random.uniform() de python et non à randrange.

Pour la recoder, comme pour beaucoup d'autres fonctions aléatoires, il vaut mieux avoir d'abord l'équivalent de random.random, qui elle est relativement simple à coder en C (mais note que la granularité va dépendre de RAND_MAX, et la qualité va dépendre de l'implémentation de rand()) :

1
2
3
double random() {
    return (double) rand() / RAND_MAX;
}

Ensuite, tu peux définir uniform comme suit (correspond à la définition de python) :

1
2
3
double uniform(double a, double b) {
    return a + (b-a) * random();
}

EDIT:

@gabbro: Il me semble que la qualité ne se dégrade pas, mais l'aléatoire devient en quelque sorte cyclique sur de très longues chaines. Typiquement, ça ne dérange en général pas vraiment pour du Monte-Carlo.

Édité par yoch

+0 -0

Oui, quasiment, je n'avais pas vu. Il te manque toutefois le cast en double (très important, et vaut mieux le faire explicitement, je ne saurais dire si le tien est transtypé au bon moment dans ton code), et je ne comprend pas bien pourquoi tu divise par RAND_MAX + 1 ? EDIT: je viens de voir que le + 1 vient de la page de comp.lang.c que tu met en lien, or c'est censé être des experts là bas. Du coup, as-tu regardé cette page ?

Édité par yoch

+0 -0
Staff

Cette réponse a aidé l'auteur du sujet

@gabbro: Il me semble que la qualité ne se dégrade pas, mais l'aléatoire devient en quelque sorte cyclique sur de très longues chaines.

Ne me demandes pas de détail, je ne suis pas spécialiste du domaine, je me suis juste pris le truc dans le nez quand j'ai fait un Monte-Carlo. :)

Typiquement, ça ne dérange en général pas vraiment pour du Monte-Carlo.

Si. On risque de laisser des états de côté. Je n'ai plus la source sous la main, mais lorsque j'ai fait du Monte-Carlo, il fallait choisir un générateur de nombre aléatoire qui soit uniforme, décorrélé et avec une longue période. De mémoire, le troisième point au moins n'est pas respecté par la fonction rand du C.

On recommande habituellement Mersenne Twister ou Well Rng. Les deux ont des implémentations libre dans à peu près tout les langages, donc pas de soucis.

Hier, dans le parc, j'ai vu une petite vieille entourée de dinosaures aviens. Je donne pas cher de sa peau.

+0 -0

Cette réponse a aidé l'auteur du sujet

La conversion en double est faite automatiquement avec la somme RAND_MAX + 1. (1. étant une constante double). Pour moi la différence entre ajouter 1 ou pas va donner un intervalle différent de valeurs mais toujours avec une granularité égale à RAND_MAX

1
2
rand()/(double)RAND_MAX; /* Intervalle [0, 1] */
rand()/(RAND_MAX+1.); /* Intervalle [0, 1( */

Sinon le code source de Mersenne/Twister en C est disponible ici, il y a une fonction qui te permet de générer un nombre aléatoire de type double sur 53 bits en combinant 2 appels à la fonction de base. Cela va te donner une granularité de 2^53.

Édité par fromvega

+0 -0

Cette réponse a aidé l'auteur du sujet

Salut,

Pour comprendre pourquoi la version modulo est biaisée, le mieux est encore de prendre un exemple. ON va admettre pour cet exemple que RAND_MAX vaut 16, et que tu veux un nombre aléatoire entre 1 et 10.

rand() peut retourner une valeur comprise entre 0 et 15, et donc un appel à randrange(1, 11) peut retourner 16 résultats différents :

  • 0 => 1
  • 1 => 2
  • 2 => 3
  • 3 => 4
  • 4 => 5
  • 5 => 6
  • 6 => 7
  • 7 => 8
  • 8 => 9
  • 9 => 10
  • 10 => 1
  • 11 => 2
  • 12 => 3
  • 13 => 4
  • 14 => 5
  • 15 => 6

Qu'est-ce qu'on remarque ? Le cycle est coupé au milieu, les valeurs 1-6 ont deux fois plus de chances de sortir que les valeurs 7-10.

Si on admet que rand() retourne effectivement un nombre aléatoire équiprobable compris entre 0 et RAND_MAX, ta fonction ne fait un tirage correct que si RAND_MAX est un multiple de l'intervalle, i.e. si (max-min)%RAND_MAX==0.

Maintenant faisons la même chose avec ta deuxième fonction :

  • 0 => 1+(0/17)*10 = 1
  • 1 => 1+(1/17)*10 = 27/17 = 1.59
  • 2 => 1+(2/17)*10 = 37/17 = 2.18
  • 3 => 47/17 = 2.76
  • 4 => 57/17 = 3.35
  • 5 => 67/17 = 3.91
  • 6 => 4.53
  • 7 => 5.12
  • 8 => 5.71
  • 9 => 6.29
  • 10 => 6.88
  • 11 => 7.47
  • 12 => 8.06
  • 13 => 8.65
  • 14 => 9.24
  • 15 => 9.82

On voit que la répartition dans l'intervalle continu 1-10 n'est pas trop mauvaise.

Si on retire le +1, on peut obtenir comme tirage maximum 15 => 1+(15/16)*10 = 166/16 = 10.38. Si en appelant randmax(1, 11) tu souhaitais un nombre entre 1 et 10, tu peux donc tomber hors bornes.

En fait, ça dépend de la définition mathématique que tu souhaites pour ta fonction. Si tu veux obtenir des nombres aléatoires réels, le +1n'a pas vraiment de sens d'après moi, car quand on écrit randrange(1, 11) on souhaite habituellement tomber sur un nombre dans l'intervalle continu semi-ouvert 1-11, i.e. on peut tomber sur 1 et sur 10.9999, mais jamais sur 11. Quand dans la plupart des langages de programmation on a un Math.random qui retourne une valeur entre 0 et 1, c'est la définition qui prévaut habituellement: on ne peut jamais tomber sur 1.

Par contre si ton but est d'obtenir des nombres aléatoires entiers, le +1 peut changer la donne. Quand on écrit randrange(1, 11), veut-on que 11 puisse faire partie des solutions possibles ou non ? Les deux sont intéressants suivant le cas: quand on jette un dé, c'est très naturel d'écrire randrange(1, 6) et que 6 fasse partie des réponses possibles; quand on veut tirer une case aléatoire dans un tableau, c'est plus simple d'écrire randrange(0, length) tout en étant certain que length ne tombera jamais; on est quitte d'oublier l'éternel -1.

Tu remarqueras que, même pour ta deuxième fonction, si on ramète ces sorties à des nombres entiers, on a le même problème, certaines valeurs ont deux fois plus de chance de tomber que d'autres. Dans la FAQ, on rappelle que l'intervalle (max-min) doit être largement inférieur à RAND_MAX pour que ce biais devienne suffisament négligeable. J'ai fait exprès pour l'exemple de prendre des petites valeurs pour que ce problème paraisse évident.

Au-delà des possibles écueils évoqués jusqu'ici, on accuse souvent rand() de ne pas être un très bon générateur de nombres aléatoires dans la plupart des implémentations de la bibliothèque C standard, pour plusieurs raisons :

  • RAND_MAX vaut souvent 65536 ou 32768 seulement, ce qui pose le problème de la condition (max-min) << RAND_MAX qui ne peut pas être terriblement remplie, et ça peut déjà se ressentir pour des tirages entre 1 et 1000.
  • rand est très souvent implémenté comme un générateur aléatoire à congruances linéaires, ce qui n'est pas un problème en soi, mais avec des paramètres particulièrement mal choisis; et comme RAND_MAX est petit, la période est forcément petite et d'autant plus avec des mauvais paramètres; donc le générateur est extrêmement prévisible.
  • Non seulement la qualité des nombres générés entre 0 et RAND_MAX n'est souvent pas terrible, mais la réalité est encore pire si on analyse les sorties bit à bit.
  • Même l'initialisation du générateur avec srand n'est souvent pas rapportée comme très bonne, car se basant uniquement sur l'heure du PC à la seconde près et avec un algorithme aussi naïf que valeur = time(NULL)%RAND_MAX. Donc même sur plusieurs exécutions successives, ça peut rester extrêmement prévisible.

Par voie de conséquence, si on veut un générateur aléatoire un minimum décent en C, on essaiera généralement de prendre une implémentation différente de celle de la bibliothèque standard. Souvent on prend un générateur Mercene Twister a.k.a. MT19937, car il est assez rapide, ne nécéssite pas trop de mémoire, et est suffisament imprévisible pour les utilisations courantes (sauf pour la cryptographie).

Pour info, un générateur MT19937 est inclus dans la bibliothèque standard du C++11.

Rajout: J'ai été multi-grillé; mais c'est pas grave.

Ma plateforme avec 23 jeux de société classiques en 6 langues et 13000 joueurs: http://qcsalon.net/ | Apprenez à faire des sites web accessibles http://www.openweb.eu.org/

+6 -0

@Renault: de fait, et à priori, je resterai sur de l'UNIX tout du long. Néanmoins, on ne sais jamais (question de rapidité, /dev/urandom est pas "lent" vu qu'il y a une lecture à faire ?).

pierre_24

Je ne connais pas bien les différents algorithmes de génération de nombres aléatoires, mais /dev/{u}random est un pseudo-fichier. Il te permet d'accéder au kernel globalement. Donc il n'y a pas de lecture disque-dur à proprement parler. Plus d'infos avec man urandom.

Le hasard n'est que le nom donné à notre ignorance et n'existerait pas pour un être ominscient., Émile Borel

+0 -0

@QuentinC : Merci pour ces explications.

Toutefois, ce que tu avance pour le +1 ne me satisfait pas. En effet, RAND_MAX fait partie de l'intervalle (dans ton cas par exemple il vaudra 15). Et on ne va donc jamais dépasser la borne (mais l'atteindre si).

Petit test :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
>>> RAND_MAX = 15
>>> for i in range(RAND_MAX + 1):
    print(i, '=>', 1 + (i / RAND_MAX) * (10-1))


0 => 1.0
1 => 1.6
2 => 2.2
3 => 2.8
4 => 3.4
5 => 4.0
6 => 4.6
7 => 5.2
8 => 5.8
9 => 6.3999999999999995
10 => 7.0
11 => 7.6
12 => 8.2
13 => 8.8
14 => 9.4
15 => 10.0

Je retiens donc plutôt l'explication que l'on veut éviter d'avoir la borne supérieure dans les résultats possible, ce qui n'est pas le cas pour random.uniform.


En ce qui concerne /dev/urandom, il me semble qu'il est très déconseillé de l'utiliser de façon massive, et je suppose que c'est lent de toutes façons.

Édité par yoch

+0 -0

Cette réponse a aidé l'auteur du sujet

@Renault: de fait, et à priori, je resterai sur de l'UNIX tout du long. Néanmoins, on ne sais jamais (question de rapidité, /dev/urandom est pas "lent" vu qu'il y a une lecture à faire ?).

pierre_24

C'est /dev/random si tu veux une certaine garantie dans la valeur de l'aléatoire produite. Et non, c'est un pseudo fichier, si ton entropie est suffisante, tu ne perdras pas plus de temps (ou à peine) à lire ces valeur qu'en les générant avec une fonction type rand().

Mais là au moins tu auras un système d'aléatoire plutôt fiable, surtout si ta machine dispose de composants pour générer de l'entropie.

Amateur de Logiciel Libre et de la distribution GNU/Linux Fedora.

+0 -0

1
2
3
>> RAND_MAX = 15
>> for i in range(RAND_MAX + 1):
    print(i, '=>', 1 + (i / RAND_MAX) * (10-1)) 

Tu ne démontres pas la même chose que moi. Soit tu fixes RAND_MAX=16 et tu boucles jusqu'à RAND_MAX et non RAND_MAX+1, soit tu divises par RAND_MAX+1 et non RAND_MAX. IL faut choisir. C'est évident que si on prend x=16 dans le calcul 1+(x/16)*10, on tombe pile sur 11.

Si RAND_MAX vaut usuellement une puissance de 2 comme 65536 ou 32768, c'est justement parce qu'elle n'est jamais atteinte par l'algorithme de congruence linéaire. Pour rappel l'algorithme de base en gros résumé c'est nombre_suivant = (nombre_actuel * A + B) % RAND_MAX avec A et B des constantes définies par l'implémentation. Sauf erreur, A et B sont censés être proche de sqrt(RAND_MAX) et premiers entre eux et avec RAND_MAX si on veut un résultat optimal: Je n'ai pas le code de la libc mais on dit partout que ces constantes A et B n'ont justement pas été bien choisis.

Édité par QuentinC

Ma plateforme avec 23 jeux de société classiques en 6 langues et 13000 joueurs: http://qcsalon.net/ | Apprenez à faire des sites web accessibles http://www.openweb.eu.org/

+0 -0
Auteur du sujet

Merci à tous, c'est très intéressant. Si je résume, ma fonction pour avoir un range est "correcte", mais le générateur de nombre aléatoire est à ch*** :p

Je vais probablement emprunter le code de Mersenne/Twister :) (merci fromvega)

@Renault:

C'est /dev/random si tu veux une certaine garantie dans la valeur de l'aléatoire produite. Et non, c'est un pseudo fichier, si ton entropie est suffisante, tu ne perdras pas plus de temps (ou à peine) à lire ces valeur qu'en les générant avec une fonction type rand().

Renault

Yep, mais la page wikipédia précise que sous Linux, l'appel à /dev/random est bloquant. À priori, j'ai besoin de générer une quantité importante de nombre aléatoire à la seconde (bon, pas des millions, mais une centaine dans le meilleur des cas), et j'ai peur que la sortie ne suive pas1. D'autant qu'à terme, c'est pour faire tourner sur un super-calculateur et plus sur ma machine, donc je ne peut pas prévoir le temps de réaction. Mais c'est une alternative.


  1. en tout cas, cat /dev/random renvoie très peu de valeur à la seconde, peut être trois caractères sur ma machine. 

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+0 -0

@Yoch: en fait, je me suis trompé sur un point. Je viens de tester rapidement (j'aurais dû le faire avant) et chez moi avec GCC 4.7 windows, RAND_MAX = 32767, et non pas 32768 comme je l'avais cru. Donc ta démonstration est correcte mais ja'vais quand même vu juste, le +1 sert bien à s'assurer qu'on ne tombera jamais sur la borne supérieure; sauf qu'il s'avère indispensable du coup. ET hop on a l'explication définitive sur pourquoi il est là.

Ma plateforme avec 23 jeux de société classiques en 6 langues et 13000 joueurs: http://qcsalon.net/ | Apprenez à faire des sites web accessibles http://www.openweb.eu.org/

+0 -0

Si RAND_MAX vaut usuellement une puissance de 2 comme 65536 ou 32768, c'est justement parce qu'elle n'est jamais atteinte par l'algorithme de congruence linéaire. Pour rappel l'algorithme de base en gros résumé c'est nombre_suivant = (nombre_actuel * A + B) % RAND_MAX avec A et B des constantes définies par l'implémentation.

QuentinC

Tu fais erreur sur ce point. Comme déjà dit, RAND_MAX fait bien partie de l'intervalle, et sa valeur est usuellement de la forme $2^n -1$ (l'implémentation étant alors nombre_suivant = (nombre_actuel * A + B) & RAND_MAX).

The rand() function shall compute a sequence of pseudo-random integers in the range [0, {RAND_MAX}] with a period of at least 232.

(The value of the {RAND_MAX} macro shall be at least 32767.)

documentation

EDIT: bon bah grilled, tant pis je laisse

Édité par yoch

+0 -0

C'est forcément en C que tu dois écrire ton code ? En C++ Mersenne-Twister est standard depuis C++11, en Fortran tu as random_number avec un vrai générateur de nombrs alératoires dans gfortran et ifort. Et avec rust, c'est dans une lib sympa !

Si tu n'as pas le choix, j'ai collecté quelques codes Monte-Carlo et j'en ai au moins un en C, avec une implémentation de générateur de nombres aléatoires. je dois pouvoir t'en passer une partie, à vérifier avec la licence.

Édité par Luthaf

Mon Github — Tuto Homebrew — Article Julia

+2 -0
Auteur du sujet

C'est forcément en C que tu dois écrire ton code ? En C++ Mersenne-Twister est standard depuis C++11, en Fortran tu as random_number avec un vrai générateur de nombrs alératoires dans gfortran et ifort. Et avec rust, c'est dans une lib sympa !

Je sais pas coder en C++ (je fais du C with class, donc autant faire du C), j'aime pas le Fortran (oui, je sais, tout les codes scientifiques de l'univers sont écrit en Fortran1, je lis le Fortran et j'aime pas ça) et je connais pas le Rust. Et je savais très bien que quelqu'un allais encore me sortir le C++11 et son (ses?) super-générateur-de-nombre-aléatoire-de-la-mort-qui-tue-disponible-en-25-milliards-de-possibilités-que-tu-sais-pas-quoi-choisir-et-les-défenseurs-de-C++-me-tapent-sur-le-système-avec-ça.2

C'est le genre de commentaires qui m'ennuient profondément. Je fais du C si je veux, et c'est pas la question. Sinon, j'aurais pas posé la question et j'aurais fait du Python (pas pour rien que je cherche l'équivalent C de la fonction Python correspondante).

Edit @Kass'Peuk: Je tape parce que on m'as encore sorti le classique "change de langage, le mien est mieux, t'as vu". Le reste est de la mauvaise foi tout à fait assumée.


  1. et j'exagère à peine. 

  2. sérieusement, sans être mathématicien et savoir manifestement quoi correspond à quoi, comment tu fais un choix correct ? Et je suis chimiste, pas mathématicien, mon cours de stat' c'est arrêté bien avant ça ! 

Édité par pierre_24

Doctorant et assistant en chimie à l'Université de NamurEx-dev' pour ZdS (a aidé à réaliser la ZEP-12 !) • Carniste cis (y parait que c'est une injure)

+1 -1
Vous devez être connecté pour pouvoir poster un message.
Connexion

Pas encore inscrit ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte