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

Я пытаюсь реализовать генератор случайных чисел с гауссовым распределением в интервале [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-м томе Кнута в TAOCP, 3-е издание, стр. 122.

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

Источник

2 ответа

Решение

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

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

Источник

Измените свое заявление 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, но я думаю, что все в порядке: по состоянию на 3.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/l2approx.pdf

Источник
Другие вопросы по тегам