Un aléatoire correct en c++

Bonjour,

Je tombe actuellement sur un problème lié à un mauvais aléatoire. Pour faire simple, j’ai une marche brownienne auto-évitante, mais elle part toujours dans la même direction.

Jusqu’ici, j’initialisais une graine et un Mersenne et des distributions,

1
2
3
4
5
6
7

mt19937_64 rng(graine);
uniform_real_distribution<real> dist_pos(-1.0l, 1.0l);
uniform_real_distribution<real> dist_p(0.l, 1.0l);

const real x_rand {dist_pos(rng)},\
y_rand {dist_pos(rng)},\
z_rand {dist_pos(rng)};

La graine est donnée à la main (c’est un entier). En pratique, ce n’est pas bon. Ce que je ne comprends vraiment pas, c’est que ça part toujours en x et y négatif, et z nul (en moyenne) – j’aurai très vite vu une erreur si ça avait été x, y et z avec le même comportement, mais là… De plus, j’ai à peu près un hasard (au sens que je n’ai pas toujours les mêmes résultats), mais avec une direction privilégiée (probablement accentué par l’auto-évitement, mais elle existe et me gêne).

Si vous avez une idée, je suis preneur…

Édit : si je fais

1
2
3

const real z_rand {dist_pos(rng)},\
y_rand {dist_pos(rng)},\
x_rand {dist_pos(rng)};

alors ça part en z et y négatifs, et x à peu près nul. Alors que c’est la même commande. Si c’est du bon aléatoire, il ne devrait pas y avoir de problème à utiliser la même graine ??? Je ne tire même pas quelques millions d’essais…

Je dis peut-être une bêtise, mais quand j’avais eu un soucis d’aléa en C, c’est parce que je lançais mon programme 10000 fois en une seconde, et du coup, la seed utilisée était trop souvent la même. A l’époque, j’avais repiqué cet algorithme.

Ca ne devrait pas jouer ici, puisque Gabbro a l’air de dire qu’il veut initialiser la seed qu’une seule fois, d’après ça :

alors ça pars en z et y négatifs, et x à peu près nul. Alors que c’est la même commande. Di c’est du bon aléatoire, il ne devrait pas y avoir de problème à utiliser la même graine ??? Je ne tire même pas quelques millions d’essais…

Je la donne à la main dans le fichier d’entrée (question de répétabilité : je peux pouvoir lancer deux fois le même calcul).

Ce qui me chiffonne vraiment, c’est que tirer x, y, z et r, ou tirer r, x, y et z donne des résultats statistiquement différent ! Dans le second cas, je retombe sur un résultat correct à vue de nez (mais auquel je ne fais guère confiance en pratique).

Bon sang,

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

const real r_rand {dist_init(rng)},\
x_rand {dist_init(rng)},\
y_rand {dist_init(rng)},\
z_rand {dist_init(rng)};
// et

const real x_rand {dist_init(rng)},\
y_rand {dist_init(rng)},\
z_rand {dist_init(rng)},\
r_rand {dist_init(rng)};

devrait pourtant bien donner le même comportement au global !

Peut-être le problème est-il causé par la graine qui est initialisée avec un seul entier, alors que l’état interne comporte 624 entiers de 32 bits ce qui pourrait donner des tirages biaisés au départ.

As-tu essayé de faire un grand nombre de tirages aléatoires "pour rien" avant de prendre ton échantillon ? Le temps que l’état interne se stabilise pour donner des tirages équiprobables.

Peut-être le problème est-il causé par la graine qui est initialisée avec un seul entier, alors que l’état interne comporte 624 entiers de 32 bits ce qui pourrait donner des tirages biaisés au départ.

J’ai modifié mon entrée pour pouvoir lui donner tout plein de trucs, mais ça ne change rien.

As-tu essayé de faire un grand nombre de tirages aléatoires "pour rien" avant de prendre ton échantillon ? Le temps que l’état interne se stabilise pour donner des tirages équiprobables.

Maintenant, oui. Et 0 tirage pour rien, 10000 ou 10 millions donnent le même résultat.

real est un typedef pour les float/double/long double (on a parfois envie de toucher à ce genre de chose sans tout réécrire dans mon domaine).

Que se passe t’il si tu remplaces mt19937_64 (réputé difficile à initialiser) par un autre générateur ?

Merci beaucoup du lien. je craignais un peu que ce soit un truc comme ça, mais je n’avais pas la source. Malheureusement, il y est écrit :

Incidentally, similar problems also affect std::minstd_rand, which uses a single 32-bit

Et en pratique, j’ai des problèmes similaire avec les autres. Je vais passer à Well RNG qui n’est pas dans la norme mais que j’ai déjà utilisé pour du Monte-Carlo et voir ce que ça donne…

Bon, l’utilisation d’un autre générateur (well rng) de nombre aléatoire marche. Le problème venait bien à priori du Mersenne, surement de l’initialisation…

Je trouve l’explication fort peu satisfaisante, puisque je faisais ce qui était proposé dans la doc, mais je vais devoir m’en contenter…

real est un typedef pour les float/double/long double (on a parfois envie de toucher à ce genre de chose sans tout réécrire dans mon domaine).

J’ai bien compris, mais je demandais à savoir pour quel type effectif le problème que tu décris arrive ? Tous ?

Tu peux essayer également PCG qui a de bonnes propriétés et est très simple a utiliser.

Mais le problème que tu décris est étrange, tu as un code minimal pour le reproduire ?

J’ai bien compris, mais je demandais à savoir pour quel type effectif le problème que tu décris arrive ? Tous ?

long double. Je viens d’essayé, et oui, j’ai le même problème en double et float.

Mais le problème que tu décris est étrange, tu as un code minimal pour le reproduire ?

Je n’ai pas réussi à en faire, il faut tout mon code de 500 lignes écrit en c++ dégueulasse (je ne suis pas informaticien) peu commenté.

Pour dire, afin de changer le générateur de nombre aléatoire, j’ai changer le code, et créé une fonction toute bête :

1
2
3
4
5
6
7
8
9

real MaClasse::alea(real mini, real maxi) {
    uniform_real_distribution<real> dist(mini, maxi);
    return dist(rng);
}

Particle Generator::new_particle_from(Particle particle, gen rng, normal_distribution<real> radius) {
    uniform_real_distribution<real> dist(-1., 1.);
    const real x_rand {alea(-1., 1.)}, y_rand {alea(-1., 1.)}, z_rand {alea(-1., 1.)}, r_rand {radius(rng)};
...

, et ça donne un résultat différent de

1
2
3
4

Particle Generator::new_particle_from(Particle particle, gen rng, normal_distribution<real> radius) {
    uniform_real_distribution<real> dist(-1., 1.);
    const real x_rand {dist(rng)}, y_rand {dist(rng)}, z_rand {dist(rng)}, r_rand {radius(rng)};
...

. Sachant que créer 3 distributions (une pour x, une pour y et une pour z) pour ne pas appeler la même ne change rien. L’ordre des appels ou l’empaquetage dans une fonction change le résultat, alors trouver un cas clairement faux et minimal, c’est compliqué. Pour tout dire, j’ai eu du pot de m’en être rendu compte…

C’est normalement que ton résultat soit différent. Dans le premier cas, tu recrées une distribution à chaque fois, dans le second tu gardes le même objet distribution.

Par « différent », il faut comprendre « statistiquement et qualitativement différent » (dans un cas, j’ai un biais dans la direction (-1, -1, 0), pour toutes graines, pas dans l’autre – j’ai surement un biais dans un coin, mais il n’apparait pas aussi clairement).

Ensuite, recréer une distribution à chaque fois (une pourx, une pour y, une pour z) dans new_particle_from ou le faire dans une fonction annexe donne des résultats (statistiquement) différent.

Dans tous les cas, je ne comprend pas pourquoi tu crée une nouvelle distribution pour chaque particule, ça semble louche puisque tu utilise visiblement la même distribution. Mais je ne vois pas trop quelle influence cela pourrait avoir…

Parce qu’avant, je ne le faisais pas, et que ça marchait pas.

Face à ce genre de bug, j’ai tendance à y aller à tâtons… Mais peut-être n’est-ce que un problème d’initialisation de Mersenne pas trivial à faire (et sans réponse facile de ce que j’ai lu).

Tu passes ton générateur par copie à new_particle_from(), donc à chaque appel ça devrait renvoyer des valeurs identiques. Passe le par référence plutôt.

Je vois que alea() n’a pas de paramètre rng, tu initialises un champ de ta classe via le paramètre rng de new_particle_from() ?

Pour alea, j’ai foutu rng dans la classe.

Pour les pointeurs, j’ai une version avec rng dans la classe (et donc pris en compte correctement lors de chaque appel, car jamais passé), et ça fait la même chose.