Генерация псевдослучайных чисел

Должен признаться, что меня никогда толком не волновли механизмы генерации случайных чисел. Ну есть себе rand() или Random, и пускай живет. Но в последнее время, особенно когда внезапно захотелось генерировать числа с нормальным распределением, пришлось немного вникнуть в тему. Поэтому этот пост – про случайные числа.

Равномерное Распределение

Когда вы вызываете какой-нибудь rand() – вы получаете псевдослучайное число с равномерным распределением. Давайте рассмотрим, почему псевдо, и почему с равномерным.

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

Кому-то может показаться, что псевдослучайно – значит плохо. Это вовсе не так. Например, если у вас есть генератор случайных картинок, вы можете сгенерировать много картинок с небольшим разрешением, а потом воссоздать понравившуюся вам картинку с большим разрешением. Главное запомнить seed value каждой картинки что, собственно, не сложно.

Теперь насчет равномерности распределения. Вообще, идея достаточно простая – когда вы генерируете число от 0 до 1, идея в том что вероятность любого числа – одинакова. Иначе говоря, если $\theta_1 < \theta_2$ , то случайная величина $Y$ имеет непрерывное равномерное вероятностное распределение на интервале $(\theta_1, \theta_2)$ если фунция плотности $Y$ выглядит вот так:

$\displaystyle f(x)=\frac{1}{\theta_{1}-\theta_{2}}, \theta_{1} \leq y \leq \theta_{2}$

Вообщем-то, самая примитивная реализация этой функции это нечто вроде $x_n = ax_{n-1}+b \ (\mathrm{mod} \ m)$ , но на практике используются чуть более сложные варианты которые дают более «качественные» случайные числа. И кому-то из вас наверное хочется закричать что «это все булшит» и «какая разница?», но когда вы делаете, например, Монте-Карло симуляции в которых случайных чисел нужно реально много, важно иметь хорошие случайные числа чтобы ваша модель ближе отражала полный хаос который творится во вселенной. (Правильное слово энтропия, кстати.)

Есть целый ряд статистических тестов для проверки качества случайных чисел. Мой подход к этому вопросу – брать лучшее, что можно получить в тех или иных случаях, например рутину ran4 из Numeical Recipes in C. Хотя, по правде сказать, даже дефолтная реализация случайных чисел в .NET – сделана «по Кнуту», и является вполне юзабельной.

Что еще сказать про равномерно распределенные случайные числа? Ну, то что среднее значение – это среднее всего диапазона, т.е.

$\displaystyle E[Y]=\frac{\theta_1+\theta_2}{2}$

А дисперсия равна

$\displaystyle V[Y]=\frac{(\theta_2-\theta_1)^2}{12}$

Последний результат особенно интересен. Попробуйте угадать, что будет если сложить 12 случайных чисел и отнять 6!

Нормальное Распределение

Нормальное распределение само по себе выглядит чуть иначе и имеет функцию плотности

$\displaystyle f(y)=\frac{1}{\sigma\sqrt{2\pi}}\exp(\frac{-(y-\mu)^2}{2\sigma^2})$

Соответственно, это распределение параметризовано средним и дисперсией, т.е. $E[Y]=\mu$ и $V[Y]=\sigma^2$ .

Грубым, но действительно работающим приближением к получению случайного нормально распределенного значения является подход, в котором мы берем, скажем, 12 случайных равномерно распределенных значений и вычитаем 6. Почему? Это возможно благодаря центральной предельной теореме, котороая говорит нам, что сумма большого количества независимых случайных величин (мы это называем i.i.d. – independent and identically distributed) имеет распределение, близкое к нормальному. Соответственно, количество случайных чисел (12) выбрано не случайно, а так чтобы «погасить» дисперсию.

Случай выше подходит разве что для простых вычислений в Excel, где нет встроенной функции RAND() с нормальным распределением. Если же говорить о программировании, то тут скорее полезно преобразование Бокса-Мюллера. Идея в том, что мы берем два случайных числа $v_1$ и $v_2$ и проверяем что они находятся в единичном круге (т.е. $v_1^2+v_2^2<1$ ). Если они подходят, то берем и возвращаем одно из этих двух чисел, наприемер $v_1$ , но не само по себе а как

$\displaystyle r = v_1^2+v_2^2 \\ f = \sqrt{-2 \frac{\log{r}}{r}} \\ \mathrm{return} \ v_1 \times f$

Доказательство этого механизма достаточно сложное, так что я его пожалуй оставлю – всегда можно почитать Википедию или Numerical Recipes. Мораль в том, что этот метод позволяет генерить большое кол-во нормально распределенных величин с определенным средним и дисперсией.

В заключение, скажу что в нормальном распределении интересно не только получение случайных величин, но также подсчет функции плотности. Интересен этот рассчет тем, что интеграл, который там фигурирует, не берется аналитически, а должен быть рассчитан численно. Соответственно, для этого есть ряд всяких реализаций с разной точностью, в т.ч. и специальные реализации для GPU.

Rate this:

3 responses to “Генерация псевдослучайных чисел”

Галина

September 12, 2013 at 12:51 am

Объясните пожалуйста,почему в ЦПТ из суммы случайных величин вычитаем 6.

1. Dmitri
  
  September 13, 2013 at 2:18 pm
  
  Tут цeль cгeнepиpoвaть чиcлo $\sim N(0,1)$ , cooтвeтcтвeннo в oбщeм cлучae нaм нужнo
  
  $\displaystyle \sqrt{\frac{12}{N}} \left( \left(\sum_{i=1}^{N} \mathtt{RAND()} \right) - \frac{N}{2} \right)$
  
  B дaннoм cлучae мы пpocтo выбpaли $N=12$ чтoбы упpocтить кopeнь cлeвa, и coтвeтcтвeннo дoлжны вычecть 6.
  
aleks

February 24, 2021 at 4:09 pm

У вас в самой первой формуле точно нет ошибки? почему там вычитается из меньшего большие? вроде должно быть “o2 – o1”, а не “o1 – o2”