Генератор гауссовских случайных чисел - PullRequest
Новая
       67

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

7 голосов
/ 13 марта 2011

Я пытаюсь реализовать генератор случайных чисел с гауссовым распределением в интервале [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));
}

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

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

Ответы [ 2 ]

8 голосов
/ 13 марта 2011

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

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

4 голосов
/ 31 января 2016

Измените свой оператор if на (s >= 1.0 || s == 0.0). Еще лучше, используйте break, как показано в следующем примере, для SIMD-гауссовского случайного числа, возвращающего комплексную пару (u, v). При этом используется генератор случайных чисел 1005 * *1003* Мерсена. Если вам нужно только одно действительное случайное число, верните только 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, но я думаю, что все в порядке: по состоянию на 3/3/2016: DLW

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

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

[3] Чарльз К. Гарретт: http://krisgarrett.net/papers/l2approx.pdf

...