Насколько плох метод rand()%N для генерации случайных чисел в небольшом диапазоне?
Собственно, все сказано в заголовке.
Часто утверждают, что этот метод плох, особенно для небольших значений N. Но хотелось бы понять, насколько он плох?
Можно ли им пользоваться или это категорически неприемлемо?
Ответы (2 шт):
Для проверки гипотезы о том, что для небольших значений N происходит перекос значений rand()%N набросал простенькую программу, которая применяет хи-квадрат критерий для разных значений N и разного количества генерируемых значений.
Саму программу можно найти тут, а результаты ее работы (VC++ 2019) — ниже.
Как видите, по сути оба способа — и старый rand()%N, и новый uniform_int_distribution + стандартный генератор — дают весьма похожие результаты.
Гипотеза о равномерном распределении на нарушается почти нигде, а если случайно и отвергается — то с равным успехом в обоих случаях.
Таким образом, как мне кажется, отвергать указанный способ не стоит. Но для каждого конкретного компилятора я бы проверял его этой программой, чтобы убедиться в том, что это и в самом деле так.
Values count = 10
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 3.60 ok 0.40 ok
3 5.9 2.60 ok 3.20 ok
4 7.8 2.00 ok 3.60 ok
5 9.5 3.00 ok 4.00 ok
6 11.1 3.20 ok 4.40 ok
7 12.6 6.80 ok 8.20 ok
8 14.1 9.20 ok 4.40 ok
9 15.5 15.20 ok 0.80 ok
Values count = 100
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.36 ok 0.00 ok
3 5.9 1.46 ok 1.34 ok
4 7.8 1.36 ok 0.56 ok
5 9.5 1.90 ok 3.70 ok
6 11.1 2.72 ok 2.12 ok
7 12.6 5.56 ok 2.62 ok
8 14.1 8.32 ok 4.96 ok
9 15.5 7.82 ok 7.10 ok
Values count = 1000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.48 ok 0.48 ok
3 5.9 0.60 ok 0.51 ok
4 7.8 4.47 ok 10.07 fail
5 9.5 3.02 ok 1.75 ok
6 11.1 3.92 ok 3.48 ok
7 12.6 9.99 ok 6.71 ok
8 14.1 19.12 fail 9.94 ok
9 15.5 5.91 ok 8.90 ok
Values count = 10000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.29 ok 1.64 ok
3 5.9 0.35 ok 0.19 ok
4 7.8 3.02 ok 3.10 ok
5 9.5 5.93 ok 7.44 ok
6 11.1 3.42 ok 3.60 ok
7 12.6 6.59 ok 3.97 ok
8 14.1 4.85 ok 4.31 ok
9 15.5 5.74 ok 10.13 ok
Values count = 100000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.18 ok 0.13 ok
3 5.9 0.21 ok 0.24 ok
4 7.8 2.57 ok 0.57 ok
5 9.5 0.55 ok 2.89 ok
6 11.1 4.19 ok 1.54 ok
7 12.6 2.71 ok 11.04 ok
8 14.1 15.91 fail 11.78 ok
9 15.5 5.87 ok 4.27 ok
Для полной самодостаточности код:
#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <random>
using namespace std;
double Xi[20][2] = // For \alpha = 0.10 && 0.05
{
{ 2.7, 3.8 }, // 1
{ 4.6, 5.9 },
{ 6.3, 7.8 },
{ 7.8, 9.5 },
{ 9.2, 11.1 },
{ 10.6, 12.6 },
{ 12.0, 14.1 },
{ 13.4, 15.5 },
{ 14.7, 16.9 },
{ 16.0, 18.3 },
{ 17.3, 19.7 },
{ 18.5, 21.0 },
{ 19.8, 22.4 },
{ 21.1, 23.7 },
{ 22.3, 25.0 },
{ 23.5, 26.3 },
{ 24.8, 27.6 },
{ 26.0, 28.9 },
{ 27.2, 30.1 },
{ 40.3, 43.8 } // > 30...
};
default_random_engine g(random_device{}());
void Experiment(int N, int Count = 10000)
{
uniform_int_distribution<> dis(0, N-1);
vector<int> r(N), u(N);
for(int i = 0; i < Count; ++i)
{
r[rand()%N]++;
u[dis(g)]++;
}
double rs = 0, us = 0;
for(int i = 0; i < N; ++i)
{
double d = double(r[i])/Count - 1./N;
rs += d*d;
d = double(u[i])/Count - 1./N;
us += d*d;
}
rs = rs * Count * N;
us = us * Count * N;
double xi = (N < 20) ? Xi[N-2][1] : Xi[19][1];
cout << setw(2) << N << " " << fixed << setprecision(1)
<< setw(5) << xi << " " << setprecision(2)
<< setw(7) << rs << ( rs < xi ? " ok " : " fail")
<< " "
<< setw(7) << us << ( us < xi ? " ok " : " fail")
<< endl;
}
int main(int argc, char * argv[])
{
srand(time(0));
for (int Count = 10; Count <= 100000; Count *= 10)
{
cout << "\n\nValues count = " << Count << "\n";
cout << " N xi(0.05) rand()%N uniform\n";
cout << "-----------------------------------------\n";
for(int N = 2; N < 10; ++N) Experiment(N,Count);
}
}
Стандарт языка не специфицирует алгоритм, используемый функцией std::rand() для генерации случайных чисел, поэтому в общем случае, не зная ни подлежащего алгоритма, ни решаемой задачи, я бы не стал делать поспешных выводов о качестве случайных чисел, генерируемых функцией rand().
Для примера рассмотрим пару реализаций функции rand().
glibc
Документация man 3 rand утверждает, что функция rand использует тот же генератор случайных чисел, что и random(), поэтому младшие биты генерируемых значений такие же случайные, как и старшие:
The versions of
rand()andsrand()in the Linux C Library use the same random number generator asrandom(3)andsrandom(3), so the lower-order bits should be as random as the higher-order bits.
Документация man 3 random поясняет, что в качестве алгоритма генерации используется «nonlinear additive feedback random number generator». Хотя конкретные детали не уточняются. Только сказано, что с помощью функции initstate() можно установить размер внутреннего состояния генератора на 8, 32, 64, 128, или даже 128 байт.
Чтобы понять, что конкретно делает random() можно посмотреть
- Исходный код функции (там в комментариях метод генерации называют «linear feedback shift register approach, employing trinomials (since there are fewer terms to sum up that way)».)
- Вопрос на enSO: glibc rand function implementation.
- Статью за авторством Peter Selinger: The GLIBC random number generator.
Изучение приведённых источников показывает, что документация тактично умалчивает, что если с помощью initstate() установить размер внутреннего состояния на восемь байт, то вместо хитрого генератора с хорошими младшими битами будет использоваться линейный конгруэнтный генератор с не очень хорошими младшими битами.
#include <iostream>
#include <stdlib.h>
using namespace std;
int main()
{
constexpr size_t state_size = 8;
char state[state_size] = {};
unsigned int seed = 1;
initstate(seed, state, state_size);
for (int N = 2; N <= 8; N *= 2)
{
cout << "N: " << N << " ";
for (int i = 0; i < 8; ++i)
{
for (int j = 0; j < N; ++j)
cout << rand() % N << " ";
cout << " ";
}
cout << "\n";
}
}
И симпатичные случайные последовательности, которые генерирует rand() % N для некоторых N:
N: 2 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
N: 4 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1
N: 8 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1
Предложенный @Harry в соседнем ответе хи-квадрат критерий данные последовательности успешно проходят:
Values count = 10
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.00 ok 1.60 ok
3 5.9 0.80 ok 0.80 ok
4 7.8 0.40 ok 0.40 ok
5 9.5 8.00 ok 8.00 ok
6 11.1 5.60 ok 3.20 ok
7 12.6 6.80 ok 12.40 ok
8 14.1 1.20 ok 9.20 ok
9 15.5 4.40 ok 11.60 ok
Values count = 100
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.00 ok 1.44 ok
3 5.9 0.32 ok 0.74 ok
4 7.8 0.00 ok 1.68 ok
5 9.5 5.80 ok 0.70 ok
6 11.1 5.84 ok 6.80 ok
7 12.6 2.34 ok 4.02 ok
8 14.1 0.16 ok 4.80 ok
9 15.5 24.74 fail 6.74 ok
Values count = 1000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.00 ok 0.32 ok
3 5.9 1.09 ok 1.23 ok
4 7.8 0.00 ok 6.37 ok
5 9.5 2.79 ok 1.27 ok
6 11.1 1.62 ok 2.90 ok
7 12.6 6.87 ok 2.34 ok
8 14.1 0.00 ok 8.00 ok
9 15.5 7.64 ok 9.17 ok
Values count = 10000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.00 ok 0.18 ok
3 5.9 2.19 ok 0.76 ok
4 7.8 0.00 ok 4.75 ok
5 9.5 8.32 ok 6.32 ok
6 11.1 2.56 ok 0.48 ok
7 12.6 2.83 ok 3.89 ok
8 14.1 0.00 ok 2.38 ok
9 15.5 4.97 ok 7.56 ok
Values count = 100000
N xi(0.05) rand()%N uniform
-----------------------------------------
2 3.8 0.00 ok 0.66 ok
3 5.9 2.62 ok 4.50 ok
4 7.8 0.00 ok 1.61 ok
5 9.5 15.34 fail 2.23 ok
6 11.1 1.73 ok 5.98 ok
7 12.6 2.51 ok 9.07 ok
8 14.1 0.00 ok 6.20 ok
9 15.5 8.17 ok 8.45 ok
MSVCRT
vc++ также использует линейный конгруэнтный генератор (Understanding the algorithm of Visual C++'s rand() function). Для улучшения свойств генерируемой последовательности часть битов отбрасывается. Из 32 битов получаемого значения отбрасываются 16 младших и один старший ((val >> 16) & 0x7fff). Это немного удлиняет периоды последовательностей, формируемых младшими разрядами.
Сформируем изображение высотой 256 пикселей и шириной image_width пикселей, где цвет каждого пикселя задаётся целым числом из отрезка [0; 255] по следующему алгоритму:
const int N = ...;
const int image_width = ...;
for (int row = 0; row < 256; ++row)
for (int col = 0; col < image_width; ++col)
image[row][col] = std::rand() % N * (255 / (N-1));
В результате получились следующие изображения.
N = 8; image_width = 509;
N = 4; image_width = 511;
N = 2; image_width = 512;
Для сравнения, если в приведённом коде заменить std::rand() на std::mt19937, то получаются такие изображения:
N = 8; image_width = 509;
N = 4; image_width = 511;
N = 2; image_width = 512;
