Избежать основного смещения rand() в моделировании Монте-Карло?

Я переписываю симуляцию Монте-Карло на C из Objective C для использования в dll из VBA/Excel. «Двигателем» расчета является создание случайного числа от 0 до 10001, которое сравнивается с переменной в районе 5000-7000. Это используется 4-800 раз за итерацию, и я использую 100000 итераций. Таким образом, получается около 50 000 000 генераций случайных чисел за прогон.

В то время как в Objective C тесты не показали предвзятости, у меня огромные проблемы с кодом C. Objective C — это надмножество C, поэтому 95% кода было скопировано, и его трудно испортить. Я прошел остальные много раз весь день вчера и сегодня, и я не нашел никаких проблем.

У меня осталась разница между arc4random_uniform() и rand() с использованием srand(), особенно из-за смещения в сторону меньших чисел от 0 до 10000. Проведенный мной тест согласуется с таким смещением от 0,5 до 2% к числам ниже примерно 5000. Любое другое объяснение состоит в том, что мой код избегал повторений, чего, я думаю, не происходит.

код очень прост ("spiller1evne" и "spiller2evne" — это числа от 5500 до 6500):

srand((unsigned)time(NULL));
for (j=0;j<antala;++j){
[..]
        for (i=1;i<450;i++){
            chance = (rand() % 10001);

[..]

             if (grey==1) {


                 if (chance < spiller1evnea) vinder = 1;
                 else vinder = 2;
            }
            else{
                if (chance < spiller2evnea) vinder = 2;
                else vinder = 1;
            }

Теперь мне не нужна настоящая случайность, вполне подойдет псевдослучайность. Мне нужно только, чтобы он был примерно равномерно распределен на кумулятивной основе (например, не имеет большого значения, если 5555 в два раза чаще выйдет, чем 5556. Имеет значение, если 5500-5599 на 5% более вероятно, чем 5600-5699 и если есть явное смещение 0,5-2% в сторону 0-4000, чем 6000-9999.

Во-первых, звучит ли правдоподобно, что rand() - моя проблема, и есть ли простая реализация, которая удовлетворяет мои низкие потребности?

РЕДАКТИРОВАТЬ: если мое подозрение правдоподобно, могу ли я использовать это:

http://www.azillionmonkeys.com/qed/random.html

Смогу ли я просто скопировать и вставить это в качестве замены (я пишу на C и использую Visual Studio, действительно новичок)?:

#include <stdlib.h>

#define RS_SCALE (1.0 / (1.0 + RAND_MAX))

double drand (void) {
    double d;
    do {
       d = (((rand () * RS_SCALE) + rand ()) * RS_SCALE + rand ()) * RS_SCALE;
    } while (d >= 1); /* Round off */
    return d;
}

#define irand(x) ((unsigned int) ((x) * drand ()))

EDIT2: ясно, что приведенный выше код работает без такой же предвзятости, поэтому я бы порекомендовал это всем, у кого есть такая же потребность в «середине дороги», как я описал выше. Это связано со штрафом, поскольку он вызывает rand() три раза. Поэтому я все еще ищу более быстрое решение.


c random vba dll arc4random
person blagstar    schedule 31.03.2015    source источник
comment
ну, простой ответ заключается в том, что вы не должны использовать rand() для серьезных симуляций Монте-Карло. rand() основан на линейной конгруэнтности, что довольно ужасно, вам следует проверить, например, другие RNG Mersenne Twister   -  person Jeanno    schedule 31.03.2015

Ответы (1)


arrow_upward
0
arrow_downward

Функция rand() генерирует int в диапазоне [0, RAND_MAX]. Если вы преобразуете это в другой диапазон с помощью оператора модуля (%), как это делает ваш исходный код, то это приведет к неравномерности, если только размер вашего целевого диапазона не будет равномерно делить RAND_MAX + 1. Это похоже на то, что вы видите.

У вас есть несколько вариантов, но если вы хотите придерживаться чего-то, основанного на rand(), я предлагаю этот вариант вашего исходного подхода:

/*
 * Returns a pseudo-random int selected from the uniform distribution
 * over the half-open interval [0, limit), provided that limit does not
 * exceed RAND_MAX.
 */
int range_rand(int limit) {
    int rand_bound = (RAND_MAX / limit) * limit;
    int r;
    while ((r = rand()) >= rand_bound) { /* empty */ }
    return r % limit;
}

Хотя в принципе количество вызовов rand(), которое генерирует каждый вызов этой функции, не ограничено, на практике среднее количество вызовов лишь немного превышает 1 для относительно небольших значений limit, а среднее значение меньше 2 для каждого значения limit. Он устраняет описанную ранее неравномерность, выбирая начальное случайное число из подмножества [0, RAND_MAX], размер которого делится без остатка на limit.

person John Bollinger    schedule 31.03.2015