Pourquoi les gens disent-ils qu'il y a un biais modulo lorsqu'on utilise un générateur de nombres aléatoires ?

Question

J'ai souvent vu cette question posée, mais je n'ai jamais vu de réponse concrète à cette question. Je vais donc en poster une ici qui, je l'espère, aidera les gens à comprendre pourquoi il y a exactement un "biais modulo" lorsqu'on utilise un générateur de nombres aléatoires, par exemple rand() en C++.

Answer 1

Así que rand() est un générateur de nombres pseudo-aléatoires qui choisit un nombre naturel entre 0 et RAND_MAX qui est une constante définie dans cstdlib (voir ce article pour un aperçu général sur rand() ).

Maintenant, que se passe-t-il si vous voulez générer un nombre aléatoire entre, disons, 0 et 2 ? Pour les besoins de l'explication, disons que RAND_MAX est 10 et je décide de générer un nombre aléatoire entre 0 et 2 en appelant rand()%3 . Cependant, rand()%3 ne produit pas les nombres entre 0 et 2 avec la même probabilité !

Quand rand() renvoie 0, 3, 6 ou 9, rand()%3 == 0 . Par conséquent, P(0) = 4/11

Quand rand() renvoie 1, 4, 7 ou 10, rand()%3 == 1 . Par conséquent, P(1) = 4/11

Quand rand() renvoie 2, 5 ou 8, rand()%3 == 2 . Par conséquent, P(2) = 3/11

Cela ne génère pas les nombres entre 0 et 2 avec une probabilité égale. Bien sûr, pour les petites fourchettes, ce n'est peut-être pas le plus gros problème, mais pour une fourchette plus large, cela pourrait fausser la distribution, en biaisant les petits nombres.

Alors, quand est-ce que rand()%n retourner une plage de nombres de 0 à n-1 avec une probabilité égale ? Lorsque RAND_MAX%n == n - 1 . Dans ce cas, avec notre hypothèse précédente rand() renvoie un nombre entre 0 et RAND_MAX avec une probabilité égale, les classes modulo de n seraient également distribuées de manière égale.

Alors comment résoudre ce problème ? Une méthode rudimentaire consiste à générer des nombres aléatoires jusqu'à ce que vous obteniez un nombre dans la fourchette souhaitée :

int x; 
do {
    x = rand();
} while (x >= n);

mais c'est inefficace pour les faibles valeurs de n puisque vous n'avez qu'un n/RAND_MAX chance d'obtenir une valeur dans votre fourchette, et vous devrez donc effectuer RAND_MAX/n appels à rand() en moyenne.

Une formule plus efficace consisterait à prendre une grande plage dont la longueur est divisible par n comme RAND_MAX - RAND_MAX % n Pour cela, continuez à générer des nombres aléatoires jusqu'à ce que vous en obteniez un qui se situe dans l'intervalle, puis prenez le module :

int x;

do {
    x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));

x %= n;

Pour de petites valeurs de n ce qui nécessitera rarement plus d'un appel à la fonction rand() .

---

Ouvrages cités et lectures complémentaires :

• Référence CPlusPlus

• Éternellement confus

---

Answer 2

La sélection aléatoire est un bon moyen d'éliminer le biais.

Mise à jour

Nous pourrions rendre le code rapide si nous recherchions un x dans l'intervalle divisible par n .

// Assumptions
// rand() in [0, RAND_MAX]
// n in (0, RAND_MAX]

int x; 

// Keep searching for an x in a range divisible by n 
do {
    x = rand();
} while (x >= RAND_MAX - (RAND_MAX % n)) 

x %= n;

La boucle ci-dessus devrait être très rapide, disons 1 itération en moyenne.

Answer 3

@user1413793 a raison sur le problème. Je ne vais pas en discuter davantage, sauf pour souligner un point : oui, pour de petites valeurs de n et de grandes valeurs de RAND_MAX le biais modulo peut être très faible. Mais l'utilisation d'un modèle induisant un biais signifie que vous devez tenir compte du biais à chaque fois que vous calculez un nombre aléatoire et choisir différents modèles pour différents cas. Et si vous faites le mauvais choix, les bogues qu'il introduit sont subtils et presque impossibles à tester en unité. Par rapport à l'utilisation d'un outil approprié (tel que arc4random_uniform ), c'est du travail supplémentaire, pas du travail en moins. Faire plus de travail pour obtenir une moins bonne solution est une mauvaise ingénierie, surtout quand il est facile de faire les choses correctement à chaque fois sur la plupart des plateformes.

Malheureusement, les mises en œuvre de la solution sont toutes incorrectes ou moins efficaces qu'elles ne devraient l'être. (Chaque solution comporte divers commentaires expliquant les problèmes, mais aucune des solutions n'a été corrigée pour les résoudre). Cela risque de perturber le chercheur de réponses occasionnel, c'est pourquoi je fournis ici une bonne implémentation connue.

Encore une fois, la meilleure solution est d'utiliser arc4random_uniform sur les plateformes qui le fournissent, ou une solution rangée similaire pour votre plateforme (telle que Random.nextInt sur Java). Il fera ce qu'il faut sans aucun coût de code pour vous. C'est presque toujours la bonne décision à prendre.

Si vous n'avez pas arc4random_uniform Vous pouvez alors utiliser la puissance de l'opensource pour voir exactement comment il est mis en œuvre au-dessus d'un RNG de plus grande portée ( ar4random dans ce cas, mais une approche similaire pourrait également fonctionner avec d'autres RNG).

Voici le Implémentation d'OpenBSD :

/*
 * Calculate a uniformly distributed random number less than upper_bound
 * avoiding "modulo bias".
 *
 * Uniformity is achieved by generating new random numbers until the one
 * returned is outside the range [0, 2**32 % upper_bound).  This
 * guarantees the selected random number will be inside
 * [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
 * after reduction modulo upper_bound.
 */
u_int32_t
arc4random_uniform(u_int32_t upper_bound)
{
    u_int32_t r, min;

    if (upper_bound < 2)
        return 0;

    /* 2**32 % x == (2**32 - x) % x */
    min = -upper_bound % upper_bound;

    /*
     * This could theoretically loop forever but each retry has
     * p > 0.5 (worst case, usually far better) of selecting a
     * number inside the range we need, so it should rarely need
     * to re-roll.
     */
    for (;;) {
        r = arc4random();
        if (r >= min)
            break;
    }

    return r % upper_bound;
}

Il est intéressant de noter le dernier commentaire du commit sur ce code pour ceux qui ont besoin d'implémenter des choses similaires :

Changez arc4random_uniform() pour calculer 2**32 % upper_bound comme -upper_bound % upper_bound . Simplifie le code et le rend identique même sur les architectures ILP32 et LP64, et aussi un peu plus rapide sur les architectures LP64 en utilisant un reste de 32 bits au lieu d'un reste de 64 bits. reste de 64 bits.

Signalé par Jorden Verwer sur tech@ ok deraadt ; pas d'objection de djm ou otto

L'implémentation Java est également facilement trouvable (voir le lien précédent) :

public int nextInt(int n) {
   if (n <= 0)
     throw new IllegalArgumentException("n must be positive");

   if ((n & -n) == n)  // i.e., n is a power of 2
     return (int)((n * (long)next(31)) >> 31);

   int bits, val;
   do {
       bits = next(31);
       val = bits % n;
   } while (bits - val + (n-1) < 0);
   return val;
 }

Answer 4

L'utilisation de modulo fait l'objet de deux reproches habituels.

• une est valable pour tous les générateurs. C'est plus facile à voir dans un cas limite. Si votre générateur a un RAND_MAX qui est 2 (ce qui n'est pas conforme à la norme C) et que vous voulez seulement 0 ou 1 comme valeur, l'utilisation de modulo générera 0 deux fois plus souvent (quand le générateur génère 0 et 2) qu'il ne générera 1 (quand le générateur génère 1). Notez que ceci est vrai dès que vous ne laissez pas tomber les valeurs, quelle que soit la correspondance que vous utilisez entre les valeurs du générateur et les valeurs souhaitées, l'une se produira deux fois plus souvent que l'autre.

• certains types de générateurs ont leurs bits les moins significatifs moins aléatoires que les autres, au moins pour certains de leurs paramètres, mais malheureusement ces paramètres ont d'autres caractéristiques intéressantes (comme le fait de pouvoir avoir RAND_MAX inférieur à une puissance de 2). Le problème est bien connu et depuis longtemps, l'implémentation des bibliothèques évite probablement le problème (par exemple, l'implémentation de rand() dans la norme C utilise ce type de générateur, mais sans les 16 bits les moins significatifs), mais certains aiment se plaindre à ce sujet et vous pouvez avoir de la malchance.

En utilisant quelque chose comme

int alea(int n){ 
 assert (0 < n && n <= RAND_MAX); 
 int partSize = 
      n == RAND_MAX ? 1 : 1 + (RAND_MAX-n)/(n+1); 
 int maxUsefull = partSize * n + (partSize-1); 
 int draw; 
 do { 
   draw = rand(); 
 } while (draw > maxUsefull); 
 return draw/partSize; 
}

pour générer un nombre aléatoire entre 0 et n évitera les deux problèmes (et il évite le débordement avec RAND_MAX == INT_MAX)

BTW, C++11 a introduit des moyens standard pour la réduction et le générateur autre que rand().

Answer 5

Comme le réponse acceptée indique, le "biais modulo" trouve ses racines dans la faible valeur de RAND_MAX . Il utilise une valeur extrêmement faible de RAND_MAX (10) pour montrer que si RAND_MAX était 10, puis que vous essayiez de générer un nombre entre 0 et 2 en utilisant %, les résultats suivants seraient obtenus :

rand() % 3   // if RAND_MAX were only 10, gives
output of rand()   |   rand()%3
0                  |   0
1                  |   1
2                  |   2
3                  |   0
4                  |   1
5                  |   2
6                  |   0
7                  |   1
8                  |   2
9                  |   0

Il y a donc 4 sorties de 0 (4/10 de chances) et seulement 3 sorties de 1 et 2 (3/10 de chances chacune).

Donc c'est biaisé. Les chiffres les plus bas ont plus de chance de sortir.

_Mais cela n'apparaît de manière évidente que lorsque RAND_MAX est petit_ . Ou plus précisément, lorsque le nombre de vos moddeurs est important par rapport au nombre de vos clients. RAND_MAX .

Une solution bien meilleure que en boucle (ce qui est incroyablement inefficace et ne devrait même pas être suggéré) est d'utiliser un PRNG avec une gamme de sortie beaucoup plus large. Le site Twister de Mersenne L'algorithme a une sortie maximale de 4 294 967 295. Ainsi, en faisant MersenneTwister::genrand_int32() % 10 à toutes fins utiles, seront distribués de manière égale et l'effet du biais modulo disparaîtra pratiquement.