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

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

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

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

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

// C #include  #include  #include  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).

Распространение: введите описание изображения здесь

Различия: введите описание изображения здесь

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

o i = (o i-3 + o i-31 ) % 2^31

или же:

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

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

x i = o i % 36 = (o i-3 + o i-31 ) % 36 = o i-3 % 36 = x i-3 ,

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

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

Мы напишем rnd i для i th вывода `random () ‘, и отметим, что:

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

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

(См. Ниже схему доказательства.)

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

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

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

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

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


random алгоритм для случая по умолчанию:

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

Для генерации rnd i : (Note: is addition mod 2 32 .)

state i = state i-31 ⊕ state i-3

rnd i = (state i - (state i 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%.

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

rnd i = (state i-31 ⊕ state i-3 - ((state i-31 + state i-3 ) mod 2)) / 2

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

= ((state i-31 - (state i-31 mod 2)) ⊕ (state i-3 - (state i-3 mod 2)) + 2) / 2

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

rnd i = rnd i-31 ⊕ rnd i-3

rnd i = rnd i-31 ⊕ rnd i-3 + 1

Как и выше, первый случай встречается в 75% случаев, предполагая, что rnd i-31 и rnd i-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) 

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

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

«[…] случайные числа не должны генерироваться с помощью метода, выбранного случайным образом. Следует использовать некоторую теорию »(DEKnuth,« Искусство компьютерного программирования », том 2)

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