Подтвердить что ты не робот

Генератор случайных чисел libc ошибочен?

Рассмотрим алгоритм проверки вероятности того, что определенное число выбрано из набора из N уникальных чисел после определенного количества попыток (например, с N = 2, какова вероятность в рулетке (без 0), которую она принимает X пытается выиграть черный?).

Правильное распределение для этого - pow (1-1/N, X-1) * (1/N).

Однако, когда я тестирую это, используя следующий код, всегда есть глубокая канава в X = 31, независимо от N, и независимо от семени.

Является ли это внутренним недостатком, который нельзя предотвратить из-за специфики реализации используемого PRNG, является ли это реальной ошибкой или я не вижу ничего очевидного?

// C

#include <sys/times.h>
#include <math.h>
#include <stdio.h>

int array[101];
void main(){

    int nsamples=10000000;
    double breakVal,diffVal;
    int i,cnt;

    // seed, but doesn't change anything
    struct tms time;
    srandom(times(&time));

    // sample
    for(i=0;i<nsamples;i++){
        cnt=1;
        do{
            if((random()%36)==0) // break if 0 is chosen
                break;
            cnt++;
        }while(cnt<100);
        array[cnt]++;
    }

    // show distribution
    for(i=1;i<100;i++){
        breakVal=array[i]/(double)nsamples; // normalize
        diffVal=breakVal-pow(1-1/36.,i-1)*1/36.; // difference to expected value
        printf("%d %.12g %.12g\n",i,breakVal,diffVal);
    }
}

Протестировано на обновленном Xubuntu 12.10 с пакетом libc6 2.15-0ubuntu20 и Intel Core i5-2500 SandyBridge, но я обнаружил это уже несколько лет назад на более старой машине Ubuntu.

Я также тестировал это на Windows 7, используя Unity3D/Mono (не уверен, какая версия Mono, хотя), и здесь канава происходит при X = 55 при использовании System.Random, в то время как Unity встроил Unity.Random не имеет видимой канавы ( по меньшей мере, не для X < 100).

Распределение: enter image description here

Различия: enter image description here

4b9b3361

Ответ 1

Это связано с тем, что функция glibc random() не является достаточно случайной. Согласно этой странице, для случайных чисел, возвращаемых random(), мы имеем:

oi = (oi-3 + oi-31) % 2^31

или

oi = (oi-3 + oi-31 + 1) % 2^31.

Теперь возьмите xi = oi % 36 и предположим, что первое уравнение выше используется (это случается с 50% шансом для каждого номера). Теперь, если xi-31=0 и xi-3!=0, тогда вероятность того, что xi=0 меньше 1/36. Это связано с тем, что 50% времени oi-31 + oi-3 будет меньше 2 ^ 31, а когда это произойдет,

xi = oi % 36 = (oi-3 + oi-31) % 36 = oi-3 % 36 = xi-3,

которая отлична от нуля. Это приводит к тому, что в канаве вы видите 31 образец после 0 отсчетов.

Ответ 2

То, что измеряется в этом эксперименте, является интервалом между успешными испытаниями эксперимента Бернулли, где успех определяется как random() mod k == 0 для некоторого k (36 в OP). К сожалению, это омрачено тем фактом, что реализация random() означает, что испытания Бернулли не являются статистически независимыми.

Мы напишем rndi для вывода ith "random()" и отметим, что:

rndi = rndi-31 + rndi-3     с вероятностью 0,75

rndi = rndi-31 + rndi-3 + 1 с вероятностью 0,25

(см. ниже рисунок для иллюстрации).

Предположим rndi-31 mod k == 0, и в настоящее время мы смотрим на rndi. Тогда это должно быть так, что rndi-3 mod k ≠ 0, потому что иначе мы бы подсчитали цикл как длину k-3.

Но (большую часть времени) (mod k): rndi = rndi-31 + rndi-3 = rndi-3 ≠ 0.

Таким образом, текущее исследование не является статистически независимым от предыдущих испытаний, а испытание 31 st после успеха намного реже будет успешным, чем в непредвзятой серии испытаний Бернулли.

Обычный совет по использованию линейно-конгруэнтных генераторов, который на самом деле не применяется к алгоритму random(), заключается в использовании битов высокого порядка вместо младших бит, поскольку бит высокого порядка "больше" случайный "(то есть меньше коррелирует с последовательными значениями). Но это не будет работать и в этом случае, потому что вышеприведенные тождества одинаково справедливы для функции high log k bits как для функции mod k == low log k bits.

Фактически, мы могли бы ожидать, что линейно-конгруэнтный генератор будет работать лучше, особенно если мы используем старшие биты вывода, потому что, хотя LCG не особенно хорош в симуляциях Монте-Карло, он не страдает от линейная обратная связь random().


random, для случая по умолчанию:

Пусть state - вектор беззнаковых длин. Инициализируйте state0...state30 с использованием семени, некоторых фиксированных значений и алгоритма микширования. Для простоты мы можем считать вектор состояния бесконечным, хотя используются только последние 31 значения, поэтому он фактически реализован как кольцевой буфер.

Для генерации rndi: (Note: добавляется mod 2 32.)

statei = statei-31 ⊕ statei-3

rndi = (statei - (statei mod 2)) / 2

Теперь обратите внимание, что:

(i + j) mod 2 = i mod 2 + j mod 2   , если i mod 2 == 0 или j mod 2 == 0

(i + j) mod 2 = i mod 2 + j mod 2 - 2, если i mod 2 == 1 и j mod 2 == 1

Если i и j равномерно распределены, первый случай будет происходить в 75% случаев, а второй случай - 25%.

Итак, путем подстановки в формулу генерации:

rndi = (statei-31 ⊕ statei-3 - ((statei-31 + statei-3) mod 2)) / 2

     = ((statei-31 - (statei-31 mod 2)) ⊕ (statei-3 - (statei-3 mod 2))) / 2 или

     = ((statei-31 - (statei-31 mod 2)) ⊕ (statei-3 - (statei-3 mod 2)) + 2) / 2

Два случая могут быть дополнительно сведены к:

rndi = rndi-31 ⊕ rndi-3

rnd i= rnd i-31 & oplus; rnd i-3 + 1

Как и выше, первый случай возникает в 75% случаев, считая, что rnd i-31 и rnd i-3 независимо выведены из равномерного распределения (которое они нет, но это разумное первое приближение).

Ответ 3

Как отмечали другие, random() не является достаточно случайным.

Использование более высоких бит вместо нижних не поможет в этом случае. Согласно руководству (man 3 rand), старые реализации rand() имели проблему в младших битах. Поэтому рекомендуется random(). Хотя в текущей реализации rand() используется тот же генератор, что и random().

Я попробовал правильное использование старого rand():

if ((int)(rand()/(RAND_MAX+1.0)*36)==0)

... и получил такую ​​же глубокую канаву при X = 31

Интересно, если я смешиваю числа rand() с другой последовательностью, я избавляюсь от канавы:

unsigned x=0;
//...

        x = (179*x + 79) % 997;
        if(((rand()+x)%36)==0)

Я использую старый Linear Congruential Generator. Я выбрал 79, 179 и 997 случайным образом из таблицы простых чисел. Это должно генерировать повторяющуюся последовательность длиной 997.

Тем не менее, этот трюк, вероятно, ввел некоторые неслучайности, некоторые след... Полученная смешанная последовательность, несомненно, не даст других статистических тестов. x никогда не принимает одинаковое значение в последовательных итерациях. Действительно, для повторения каждого значения требуется ровно 997 итераций.

'' [..] случайные числа не должны генерироваться с помощью метода, выбранного случайным образом. Следует использовать некоторую теорию "(Д.Е.Кнут," Искусство компьютерного программирования ", том 2)

Для моделирования, если вы хотите быть уверенным, используйте Mersenne Twister