C: Генерация Случайных Чисел-Что (Если Что-То) Не Так С Этим


Для простого моделирования в C мне нужно сгенерировать экспоненциальные случайные величины. Я помню, что где-то читал (но сейчас не могу найти, и не помню почему), что использование функции rand() для генерации случайных целых чисел в фиксированном диапазоне будет генерировать неравномерно распределенные целые числа. Из-за этого мне интересно, может ли этот код иметь подобную проблему:

//generate u ~ U[0,1]
u = (   (double)rand() / ((double)(RAND_MAX));
//inverse of exponential CDF to get exponential random variable
expon = -log(1-u) * mean;

Спасибо!

5   4  

5 ответов:

Проблема со случайными числами в фиксированном диапазоне заключается в том, что многие люди делают это для чисел от 100 до 200, например:

100 + rand() % 100

Это не однообразно. Но при этом он является (или, по крайней мере, достаточно близок к единому):

u = 100 + 100 * ((double)rand() / ((double)(RAND_MAX));
Так как это то, что вы делаете, вы должны быть в безопасности.

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

RAND_MAX обычно равен 32k, в то время как использование LCG rand() генерирует псевдослучайные 32-разрядные числа. Таким образом, отсутствие однородности, а также низкая периодичность, как правило, остаются незамеченными.

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

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

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

Что касается "нормально", всегда возможно, что у вас есть 64-битный генератор случайных чисел, где у двойника обычно есть 53-битная мантисса. В этом случае у вас может возникнуть та же проблема, что и с зажимом диапазона целыми числами.

Нет, ваш алгоритм будет работать; он использует функцию модуля, которая делает вещи несовершенно.
Одна проблема заключается в том, что, поскольку он квантован, время от времени он будет генерировать именно RAND_MAX, и вы будете просить log(1-1). Я бы рекомендовал по крайней мере (rand() + 0.5)/(RAND_MAX+1), если не лучший источник, такой как drand48().

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