Почему люди говорят, что есть смещение по модулю при использовании генератора случайных чисел?


Я видел, что этот вопрос задают много, но никогда не видел истинного конкретного ответа на него. Поэтому я собираюсь опубликовать здесь один, который, надеюсь, поможет людям понять, почему именно существует "смещение по модулю" при использовании генератора случайных чисел, например rand() в C++.

9 247

9 ответов:

так rand() это генератор псевдослучайных чисел, который выбирает натуральное число между 0 и RAND_MAX, которая является константой, определенной в cstdlib (см. статьи для общего обзора на rand()).

теперь что произойдет, если вы хотите создать случайное число между скажем 0 и 2? Ради объяснения, скажем RAND_MAX 10, и я решил создать случайное число между 0 и 2, позвонив rand()%3. Однако,rand()%3 не производят числа от 0 до 2 с равной вероятностью!

, когда rand() возвращает 0, 3, 6, или 9,rand()%3 == 0. Поэтому P(0) = 4/11

, когда rand() возвращает 1, 4, 7, или 10,rand()%3 == 1. Поэтому P(1) = 4/11

, когда rand() возвращает 2, 5, или 8,rand()%3 == 2. Поэтому P(2) = 3/11

это не создает числа между 0 и 2 с равными вероятность. Конечно, для небольших диапазонов это может быть не самая большая проблема, но для большего диапазона это может исказить распределение, смещая меньшие числа.

когда rand()%n вернуть диапазон чисел от 0 до n-1 с равной вероятностью? Когда RAND_MAX%n == n - 1. В этом случае, наряду с нашим предыдущим предположением rand() возвращает число между 0 и RAND_MAX С равной вероятностью классы по модулю n также будут равномерно распределены.

так как мы решим эту проблему? Грубый способ заключается в том, чтобы продолжать генерировать случайные числа, пока вы не получите число в желаемом диапазоне:

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

но это неэффективно для низких значений n, так как у вас есть только n/RAND_MAX шанс получить значение в вашем диапазоне, и поэтому вам нужно будет выполнить RAND_MAX/n звонки rand() в среднем.

более эффективным формульным подходом было бы взять некоторый большой диапазон с длиной, кратной n, как RAND_MAX - RAND_MAX % n, хранить случайные числа, пока вы не получите тот, который лежит в диапазоне, а затем взять модуль:

int x;

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

x %= n;

для малых значений n, это редко потребует более одного вызова rand().


процитированные работы и дальнейшее чтение:


продолжайте выбирать случайный-это хороший способ удалить смещение.

обновление

мы могли бы сделать код быстро, если мы ищем x в диапазоне, кратном 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;

приведенный выше цикл должен быть очень быстрым, скажем, 1 итерация в среднем.

@user1413793 правильно о проблеме. Я не собираюсь обсуждать это дальше, за исключением одного момента: да, для небольших значений n и большие значения RAND_MAX, смещение по модулю может быть очень маленьким. Но Использование шаблона, вызывающего смещение, означает, что вы должны учитывать смещение каждый раз, когда вы вычисляете случайное число и выбираете разные шаблоны для разных случаев. И если вы сделаете неправильный выбор, ошибки, которые он вводит, являются тонкими и почти невозможными для модульного тестирования. По сравнению с просто используя правильный инструмент (например,arc4random_uniform), это дополнительная работа, а не меньше работы. Делать больше работы и получать худшее решение-это ужасная инженерия, особенно когда делать это правильно каждый раз легко на большинстве платформ.

к сожалению, реализации решения все неверны или менее эффективны, чем они должны быть. (Каждое решение имеет различные комментарии, объясняющие проблемы, но ни одно из решений не было исправлено для их решения.) Это скорее всего смутит случайный искатель ответов, поэтому я предоставляю здесь известную хорошую реализацию.

опять же, лучшее решение-просто использовать arc4random_uniform на платформах, которые предоставляют его, или аналогичное ранжированное решение для вашей платформы (например,Random.nextInt на Java). Он будет делать правильные вещи без каких-либо затрат код для вас. Это почти всегда правильное решение.

если у вас нет arc4random_uniform, то вы можете использовать силу opensource, чтобы увидеть, как именно это происходит реализовано поверх более широкого диапазона ГСЧ (ar4random в этом случае, но аналогичный подход может также работать поверх других ГСЧ).

здесь реализация 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;
}

стоит отметить последний комментарий commit к этому коду для тех, кому нужно реализовать подобные вещи:

изменить arc4random_uniform() для вычисления 2**32 % upper_bound'' as -объектом upper_bound % объектом upper_bound". Упрощает код и делает его то же самое на обоих Архитектуры ILP32 и LP64, а также немного быстрее Архитектуры LP64 с использованием 32-разрядного остатка вместо 64-разрядного остаток.

указал Джорден Вервер на технологии@ ok deraadt; никаких возражений от djm или otto

реализация Java также легко найти (см. предыдущую ссылку):

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;
 }

определение

Смещение По Модулю является врожденным смещением при использовании арифметики по модулю для уменьшения выходного набора до подмножества входного набора. В общем случае смещение существует всякий раз, когда отображение между входным и выходным набором не распределено равномерно, как в случае использования арифметики по модулю, когда размер выходного набора не является делителем размера входного набора.

это смещение особенно трудно избежать в вычислениях, где числа представлены как строки битов: 0s и 1S. Найти по-настоящему случайных источников случайности тоже крайне сложно, но это выходит за рамки данного обсуждения. в течение оставшейся части этого ответа, предположим, что существует неограниченный источник действительно случайных битов.

давайте рассмотрим моделирование рулона штампа (от 0 до 5) с использованием этих случайных битов. Есть 6 возможностей, поэтому нам нужно достаточно битов, чтобы представить число 6, которое составляет 3 бита. К сожалению, 3 случайных бита дают 8 возможных результатов:

000 = 0, 001 = 1, 010 = 2, 011 = 3
100 = 4, 101 = 5, 110 = 6, 111 = 7

мы можем уменьшить размер результирующего набора ровно до 6, взяв значение по модулю 6, однако это представляет смещение по модулю: 110 дает 0, а 111 возвращает 1. этот штамп загружен.

Возможные Решения

подход 0:

вместо того, чтобы полагаться на случайные биты, в теории можно нанять небольшую армию, чтобы бросить кости все день и запишите результаты в базу данных, а затем используйте каждый результат только один раз. Это примерно так же практично, как это звучит, и, скорее всего, не даст действительно случайных результатов в любом случае (каламбур).

подход 1:

вместо использования модуля наивное, но математически правильное решение состоит в том, чтобы отбросить результаты, которые дают 110 и 111 и просто попробуйте еще раз с 3 новыми битами. К сожалению, это означает, что есть 25% шанс на каждый рулон, что потребуется повторный рулон, включая каждый из повторных рулонов сами. Это явно непрактично для всех, кроме самого тривиального использования.

подход 2:

Используйте больше битов: вместо 3 битов используйте 4. Это дает 16 возможных результатов. Конечно, повторная прокатка в любое время, когда результат больше 5, ухудшает ситуацию (10/16 = 62,5%), так что это само по себе не поможет.

обратите внимание, что 2 * 6 = 12

звучит хорошо на первый взгляд, но давайте проверим математику:

4 discarded results / 16 possibilities = 25%

в этом случае 1 дополнительный бит не помог на всех!

этот результат неудачен, но давайте попробуем еще раз с 5 битами:

32 % 6 = 2 discarded results; and
2 discarded results / 32 possibilities = 6.25%

определенное улучшение, но не достаточно хорошо во многих практических случаях. Этот хорошая новость,добавление большего количества битов никогда не увеличит шансы на необходимость отбросить и повторно свернуть. Это справедливо не только для кубиков, но и во всех случаях.

как показала однако добавление 1 дополнительного бита может ничего не изменить. на самом деле, если мы увеличим наш рулон до 6 бит, вероятность остается 6.25%.

это вызывает 2 дополнительных вопроса:

  1. если мы добавим достаточно битов, есть ли гарантия, что вероятность a отбросы уменьшатся?
  2. сколько бит достаточно в общем случае?

Общее Решение

к счастью, ответ на первый вопрос-да. Проблема с 6 заключается в том, что 2^x mod 6 переворачивается между 2 и 4, которые по совпадению кратны 2 друг от друга, так что для четного x > 1,

[2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1)

таким образом, 6 является скорее исключением, чем правилом. Можно найти большие модули, которые дают последовательные степени 2 таким же образом, но в конечном итоге это должно обернуться, и вероятность сброса будет уменьшена.

без предоставления дополнительных доказательств, в общем, используя удвоить количество бит требуется обеспечит меньший, как правило, незначительно., шанс на выброс.

доказательство концепции

вот пример программы, которая использует libcrypo OpenSSL для предоставления случайных байтов. При компиляции, быть обязательно ссылка на библиотеку с -lcrypto который большинство людей должны иметь в наличии.

#include <iostream>
#include <assert.h>
#include <limits>
#include <openssl/rand.h>

volatile uint32_t dummy;
uint64_t discardCount;

uint32_t uniformRandomUint32(uint32_t upperBound)
{
    assert(RAND_status() == 1);
    uint64_t discard = (std::numeric_limits<uint64_t>::max() - upperBound) % upperBound;
    uint64_t randomPool = RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));

    while(randomPool > (std::numeric_limits<uint64_t>::max() - discard)) {
        RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
        ++discardCount;
    }

    return randomPool % upperBound;
}

int main() {
    discardCount = 0;

    const uint32_t MODULUS = (1ul << 31)-1;
    const uint32_t ROLLS = 10000000;

    for(uint32_t i = 0; i < ROLLS; ++i) {
        dummy = uniformRandomUint32(MODULUS);
    }
    std::cout << "Discard count = " << discardCount << std::endl;
}

я призываю играть с MODULUS и ROLLS значения, чтобы увидеть, сколько ре-роллы на самом деле происходит в большинстве случаев. Скептически настроенный человек может также пожелать сохранить вычисленные значения в файл и проверить, что распределение выглядит нормально.

есть две обычные жалобы с использованием модуля.

  • одна действителен для всех генераторов. Это легче увидеть в предельном случае. Если ваш генератор имеет RAND_MAX, который равен 2 (что не соответствует стандарту C), и вы хотите только 0 или 1 в качестве значения, используя modulo будет генерировать 0 в два раза чаще (когда генератор генерирует 0 и 2), как он будет генерировать 1 (когда генератор генерирует 1). Обратите внимание, что это верно, как только вы не бросайте значения, что сопоставление, которое вы используете от значений генератора до нужного, будет происходить в два раза чаще, чем другое.

  • некоторые генераторы имеют свои менее значимые биты менее случайными, чем другие, по крайней мере для некоторых из их параметров, но, к сожалению, эти параметры имеют другую интересную характеристику (такая имеет возможность иметь RAND_MAX один меньше, чем степень 2). Проблема хорошо известна и в течение длительного времени реализации библиотеки, наверное, не проблема (например, пример реализации rand() в стандарте C использует такой генератор, но отбрасывает 16 менее значимых бит), но некоторые любят жаловаться на это, и вам может не повезти

использовать что-то вроде

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; 
}

чтобы сгенерировать случайное число между 0 и n, вы избежите обеих проблем (и это позволит избежать переполнения RAND_MAX = = INT_MAX)

кстати, C++11 ввел стандартные способы сокращения и другие генераторы, чем Ранд.)(

решение Марка (принятое решение) почти идеально.

int x;

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

x %= n;

отредактировано 25 марта '16 в 23:16

Марк Амери 39k21170211

однако он имеет оговорку, которая отбрасывает 1 допустимый набор результатов, отбрасываемых во всех сценариях, где значение RAND_MAX (RM) на 1 меньше кратного N.

т. е., когда число значений, которые будут отброшены как недопустимые (I), равно N, то они являются фактически допустимый набор, а не недопустимый набор.

например:

RM = 255
N = 4

Discard X => RM - RM % N 

When X => 252, Discarded values = 252, 253, 254, 255

Number of discarded Values (I) = RM % N + 1

как вы можете видеть в Примере количество отброшенных значений = 4, когда количество отброшенных значений = N, то набор допустим для использования.

если мы опишем разницу между значениями N и RM как D, то есть:

D = (RM - N)

затем, когда значение D становится меньше, процент ненужных повторных рулонов из-за этого метода увеличивается при каждом естественном мультипликативном. (Так когда же RAND_MAX не равно простому числу это допустимая проблема)

например:

RM=255 , N=2 Then: D = 253, Lost percentage = 0.78125%

RM=255 , N=4 Then: D = 251, Lost percentage = 1.5625%
RM=255 , N=8 Then: D = 247, Lost percentage = 3.125%
RM=255 , N=16 Then: D = 239, Lost percentage = 6.25%
RM=255 , N=32 Then: D = 223, Lost percentage = 12.5%
RM=255 , N=64 Then: D = 191, Lost percentage = 25%
RM=255 , N= 128 Then D = 127, Lost percentage = 50%

чтобы отрицать это, мы можем сделать простую поправку, как показано здесь:

 int x;

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

 x %= n;

это обеспечивает более общую версию формулы, которая учитывает дополнительные особенности использования модуля для определения ваших максимальных значений.

примеры использования небольшого значения для RAND_MAX, которое является мультипликативным от N.

Марк Ориджинал Версия:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X >= (RAND_MAX - ( RAND_MAX % n ) )
When X >= 2 the value will be discarded, even though the set is valid.

Исправленный Вариант:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X > (RAND_MAX - ( ( RAND_MAX % n  ) + 1 ) % n )
When X > 3 the value would be discarded, but this is not a vlue in the set RAND_MAX so there will be no discard.

кроме того, в случае, когда N должно быть число значений в RAND_MAX; в этом случае вы можете установить N = RAND_MAX +1, Если RAND_MAX = INT_MAX.

по циклу вы можете просто использовать N = 1, и любое значение X будет принято, однако, и поставить оператор IF для вашего конечного множителя. Но, возможно, у вас есть код, который может иметь вескую причину для возврата 1, когда функция вызывается с n = 1...

поэтому может быть лучше использовать 0, который обычно обеспечивает ошибку Div 0, когда вы хотите иметь n = RAND_MAX+1

IE:

int x;

if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

    x %= n;
} else {
    x = rand();
}

оба этих решения решают проблему с ненужным отбрасыванием допустимых результатов, которые будут возникать, когда RM+1 является продуктом n.

вторая версия также охватывает сценарий граничного случая, когда вам нужно N, чтобы равняться общему возможному набору значений, содержащихся в RAND_MAX.

в модифицированный подход в обоих случаях одинаков и позволяет получить более общее решение необходимости предоставления допустимых случайных чисел и минимизации отброшенных значений.

повторю:

основное общее решение, которое расширяет пример Марка:

 int x;

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

 x %= n;

расширенное общее решение, которое позволяет один дополнительный сценарий RAND_MAX+1 = n:

int x;

if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

    x %= n;
} else {
    x = rand();
}

С RAND_MAX стоимостью 3 (на самом деле он должен быть намного выше, но смещение все равно будет существовать) из этих расчетов имеет смысл, что есть смещение:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

в этом случае % 2 это то, что вы не должны делать, когда вы хотите случайное число между 0 и 1. Вы можете получить случайное число между 0 и 2 делать % 3 хотя бы потому, что в этом случае: RAND_MAX это несколько 3.

другой способ

там гораздо проще, но, чтобы добавить к другим ответам, вот мое решение, чтобы получить случайное число между 0 и n - 1, так что n разные возможности, без предвзятости.

  • количество битов (не байтов), необходимых для кодирования количество возможностей-это количество битов случайных данных, которые вам понадобятся
  • кодировать число из случайных битов
  • если это число >= n перезагрузите (не по модулю).

действительно случайные данные нелегко получить, так зачем использовать больше битов, чем нужно.

Ниже приведен пример в Smalltalk, используя кэш битов из генератора псевдослучайных чисел. Я не эксперт по безопасности, так что используйте на свой страх и риск.

next: n

    | bitSize r from to |
    n < 0 ifTrue: [^0 - (self next: 0 - n)].
    n = 0 ifTrue: [^nil].
    n = 1 ifTrue: [^0].
    cache isNil ifTrue: [cache := OrderedCollection new].
    cache size < (self randmax highBit) ifTrue: [
        Security.DSSRandom default next asByteArray do: [ :byte |
            (1 to: 8) do: [ :i |    cache add: (byte bitAt: i)]
        ]
    ].
    r := 0.
    bitSize := n highBit.
    to := cache size.
    from := to - bitSize + 1.
    (from to: to) do: [ :i |
        r := r bitAt: i - from + 1 put: (cache at: i)
    ].
    cache removeFrom: from to: to.
    r >= n ifTrue: [^self next: n].
    ^r

как принято отвечать указывает ,что "смещение по модулю" имеет свои корни в низком значении RAND_MAX. Он использует чрезвычайно малое значение RAND_MAX (10) чтобы показать, что если RAND_MAX было 10, то вы попытались сгенерировать число от 0 до 2 с помощью%, будут получены следующие результаты:

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

таким образом, есть 4 выхода 0-х (4/10 шанс) и только 3 выхода 1 и 2 (3/10 шансы каждый).

так что это необъективно. Меньшие числа, имеют больше шансов выходит.

но это только проявляется так очевидно, когда RAND_MAX маленький. Или, более конкретно, когда номер вашего моддинга является большим по сравнению с RAND_MAX.

гораздо лучшее решение, чем циклы (что безумно неэффективно и даже не должно быть предложено) заключается в использовании PRNG с гораздо большим диапазоном выходных данных. Элемент Мерсенн Твистер алгоритм имеет максимальный выход 4,294,967,295. Как такое делают MersenneTwister::genrand_int32() % 10 для всех намерений и целей, будут равномерно распределены и эффект смещения по модулю будет все, но исчезнет.

Я только что написал код для метода непредвзятого подбрасывания монет фон Неймана, который теоретически должен устранить любое смещение в процессе генерации случайных чисел. Дополнительную информацию можно найти по адресу (http://en.wikipedia.org/wiki/Fair_coin)

int unbiased_random_bit() {    
    int x1, x2, prev;
    prev = 2;
    x1 = rand() % 2;
    x2 = rand() % 2;

    for (;; x1 = rand() % 2, x2 = rand() % 2)
    {
        if (x1 ^ x2)      // 01 -> 1, or 10 -> 0.
        {
            return x2;        
        }
        else if (x1 & x2)
        {
            if (!prev)    // 0011
                return 1;
            else
                prev = 1; // 1111 -> continue, bias unresolved
        }
        else
        {
            if (prev == 1)// 1100
                return 0;
            else          // 0000 -> continue, bias unresolved
                prev = 0;
        }
    }
}