Генератор случайных чисел Гаусса

Я пытаюсь реализовать генератор случайных чисел с распределением по Гауссу в интервале [0,1].

float rand_gauss (void) {
  float v1,v2,s;

  do {
    v1 = 2.0 * ((float) rand()/RAND_MAX) - 1;
    v2 = 2.0 * ((float) rand()/RAND_MAX) - 1;

    s = v1*v1 + v2*v2;
  } while ( s >= 1.0 );

  if (s == 0.0)
    return 0.0;
  else
    return (v1*sqrt(-2.0 * log(s) / s));
}

Это в значительной степени прямая реализация алгоритма Кнута во 2-м томе 3-го издания TAOCP, стр. 122.

Проблема в том, что rand_gauss() иногда возвращает значения за пределами интервала [0,1].


c random gaussian taocp
person Community    schedule 13.03.2011    source источник
comment
Гауссиана неограничена. Я что-то упускаю?   -  person Dr. belisarius    schedule 13.03.2011
comment
есть дисперсия и среднее значение, я беру среднее значение за 0, а дисперсию ^ 2 за 1, то есть стандартное нормальное распределение.   -  person    schedule 13.03.2011
comment
@nvm: стандартное нормальное распределение может с некоторой вероятностью принимать любое значение от -бесконечности до бесконечности; нет ограничений по диапазону результатов.   -  person Jeremiah Willcock    schedule 13.03.2011
comment
Вы правы, я воспринял это как рецепт и особо не думал об этом. Эпический провал :|. Спасибо!   -  person    schedule 13.03.2011
comment
Кто-нибудь опубликуйте свой комментарий как ответ, чтобы я мог пометить его как решенный.   -  person    schedule 13.03.2011
comment
@nvm Теперь, когда вы просветлены, я думаю, что лучший способ - опубликовать ответ самостоятельно и принять его, возможно, упомянув комментарии.   -  person Dr. belisarius    schedule 13.03.2011
comment
Есть ли причина, по которой вы используете float, а не double? Обычно это плохая идея..   -  person R.. GitHub STOP HELPING ICE    schedule 13.03.2011
comment
@R.. нет особой причины, не хочешь уточнить? Почему это плохая идея? Мне достаточно 6-7 знаков точности.   -  person    schedule 13.03.2011
comment
Арифметика для float и double почти наверняка стоит одинаково, плюс вы все равно конвертируете туда и обратно в double, когда вызываете log и sqrt.   -  person R.. GitHub STOP HELPING ICE    schedule 13.03.2011

Ответы (2)


arrow_upward
8
arrow_downward

Кнут описывает полярный метод на стр. 122 2-го тома TAOCP. Этот алгоритм генерирует нормальное распределение со средним значением = 0 и стандартным отклонением = 1. Но вы можете настроить это, умножив желаемое стандартное отклонение и добавив желаемое среднее значение.

Возможно, вам будет интересно сравнить свой код с другой реализацией полярного метода в C-FAQ.

person Mike Sherrill 'Cat Recall'    schedule 13.03.2011

arrow_upward
4
arrow_downward

Измените оператор if на (s >= 1.0 || s == 0.0). Еще лучше использовать break, как показано в следующем примере, для SIMD-генерации гауссовых случайных чисел, возвращающих сложную пару (u, v). При этом используется генератор случайных чисел скручивания Мерсенна dsfmt(). Если вам нужно только одно реальное случайное число, верните только u и сохраните v для следующего прохода.

inline static void randn(double *u, double *v)
{
double s, x, y;    // SIMD Marsaglia polar version for complex u and v
while (1){
   x = dsfmt_genrand_close_open(&dsfmt) - 1.; 
   y = dsfmt_genrand_close_open(&dsfmt) - 1.;
   s = x*x + y*y;
   if (s < 1) break;
}
 s = sqrt(-2.0*log(s)/s);
*u = x*s; *v = y*s;
return;
}

Этот алгоритм удивительно быстр. Время выполнения для вычисления двух случайных чисел (u, v) для четырех разных генераторов гауссовских случайных чисел составляет:

Times for delivering two Gaussian numbers (u + iv)
i7-2600K @ 4GHz, gcc -Wall -Ofast -msse2 ..

gsl_ziggurat                        =       20.3 (ns) 
Box-Muller                          =       78.8 (ns) 
Box-Muller with fast_sin fast_cos   =       28.1 (ns) 
SIMD Marsaglia polar                =       35.0 (ns)

Полиномиальные процедуры fast_sin и fast_cos Чарльза К. Гарретта ускоряют вычисления Бокса-Мюллера в 2,9 раза, используя вложенную полиномиальную реализацию функций cos() и sin(). SIMD Box Muller и полярные алгоритмы, безусловно, конкурентоспособны. Также их можно легко распараллелить. Используя gcc -Ofast -S, дамп ассемблерного кода показывает, что квадратный корень — это SIMD SSE2: sqrt --> sqrtsd %xmm0, %xmm0

Комментарий: получить точные тайминги с помощью gcc5 действительно сложно и неприятно, но я думаю, что это нормально: по состоянию на 03.02.2016: DLW

[1] Ссылка по теме: возврат указателя массива c malloc в cython

[2] Сравнение алгоритмов, но не обязательно для версий SIMD: http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

[3] Чарльз К. Гарретт: http://krisgarrett.net/papers/l2приблизительно.pdf

person W9DKI    schedule 31.01.2016
comment
вот сложная часть, чтобы понять правильно: gcc -Wall -Ofast -msse2 -frename-registers -malign-double -fno-strict-aliasing -DHAVE_SSE2=1 -DDSFMT_MEXP=19937 -o randn_test dSFMT.c randn_test.c -lm - - person W9DKI; 31.01.2016
comment
Я полагаю, @W9DKI, вы по ошибке создали дублирующую учетную запись. Когда вы войдете в следующий раз, войдите в систему так же, как и в первый раз, а не каким-либо другим способом. Также отметьте для слияния учетных записей. См. этот stackoverflow.com/help/merging-accounts. - person Bhargav Rao; 01.02.2016
comment
Джонатан, спасибо за замену опечатки s=0 на s==0. Вместо использования while(1) for(;;) было бы немного сложнее, но gcc -Ofast компилирует оба варианта в одно и то же. - person W9DKI; 07.02.2016
comment
здесь явно пропущен множитель 2, т. е. x = 2.*dsfmt_genrand_close_open(&dsfmt) - 1,0; у = 2.*dsfmt_genrand_close_open(&dsfmt) - 1,0 - person W9DKI; 22.03.2016
comment
x = dsfmt_genrand_close_open(&dsfmt) - 1.; y = dsfmt_genrand_close_open(&dsfmt) - 1.; - person W9DKI; 22.03.2016