Полярный метод Марсальи

Marsaglia полярный метод является номер выборки псевдослучайного методом генерации пары независимых стандартных нормальных случайных величин. Хотя алгоритм Зикгурата превосходит преобразование Бокса – Мюллера, он даже более эффективен.

Стандартные нормальные случайные величины часто используются в информатике, вычислительной статистике и, в частности, в приложениях метода Монте-Карло.

Полярный метод работает путем выбора случайных точек ( x,  y) в квадрате −1 lt;  x  lt;1, −1 lt;  y  lt;1, пока

${\ Displaystyle 0 lt;s = х ^ {2} + y ^ {2} lt;1, \,}$

а затем возвращая требуемую пару нормальных случайных величин как

${\ displaystyle x {\ sqrt {\ frac {-2 \ ln (s)} {s}}} \,, \ \ y {\ sqrt {\ frac {-2 \ ln (s)} {s}}},}$

или, что то же самое,

${\ displaystyle {\ frac {x} {\ sqrt {s}}} {\ sqrt {-2 \ ln (s)}} \,, \ \ {\ frac {y} {\ sqrt {s}}} { \ sqrt {-2 \ ln (s)}},}$

где и представляют собой косинус и синус угла, который вектор ( x, y) образует с осью x. ${\ Displaystyle х / {\ sqrt {s}}}$${\ displaystyle y / {\ sqrt {s}}}$

• 1 Теоретическая основа
• 2 История
• 3 Практические соображения
• 4 Реализация
• 5 ссылки

Основную теорию можно резюмировать следующим образом:

Если u равномерно распределен в интервале 0 ≤  u  lt;1, то точка (cos (2π u), sin (2π u)) равномерно распределена на единичной окружности x ²  +  y ²  = 1, и умножая эту точку на независимая случайная величина ρ, распределение которой

${\ displaystyle \ Pr (\ rho lt;a) = \ int _ {0} ^ {a} re ^ {- r ^ {2} / 2} \, dr}$

поставит точку

${\ Displaystyle \ влево (\ ро \ соз (2 \ пи и), \ ро \ грех (2 \ пи и) \ справа)}$

координаты которых совместно распределены как две независимые стандартные нормальные случайные величины.

Эта идея восходит к Лапласу, которому Гаусс приписывает открытие вышеуказанного

${\ displaystyle I = \ int _ {- \ infty} ^ {\ infty} e ^ {- x ^ {2} / 2} \, dx}$

извлекая квадратный корень из

${\ displaystyle I ^ {2} = \ int _ {- \ infty} ^ {\ infty} \ int _ {- \ infty} ^ {\ infty} e ^ {- (х ^ {2} + y ^ {2 }) / 2} \, dx \, dy = \ int _ {0} ^ {2 \ pi} \ int _ {0} ^ {\ infty} re ^ {- r ^ {2} / 2} \, dr \, d \ theta.}$

Преобразование в полярные координаты показывает, что θ равномерно распределен (постоянная плотность) от 0 до 2π и что радиальное расстояние r имеет плотность

${\ displaystyle re ^ {- r ^ {2} / 2}. \,}$

( r ² имеет соответствующее распределение хи-квадрат. )

Этот метод создания пары независимых стандартных нормальных переменных путем радиального проецирования случайной точки на единичной окружности на расстояние, задаваемое квадратным корнем из переменной хи-квадрат-2, называется полярным методом генерации пары нормальных случайных величин.,

Прямое применение этой идеи,

${\ displaystyle x = {\ sqrt {-2 \ ln (u_ {1})}} \ cos (2 \ pi u_ {2}), \ quad y = {\ sqrt {-2 \ ln (u_ {1}))}} \ sin (2 \ pi u_ {2})}$

называется преобразованием Бокса – Маллера, в котором переменная хи обычно генерируется как

${\ displaystyle {\ sqrt {-2 \ ln (u_ {1})}},}$

но это преобразование требует функций логарифма, квадратного корня, синуса и косинуса. На некоторых процессорах косинус и синус одного и того же аргумента могут быть вычислены параллельно с помощью одной инструкции. В частности, для машин на базе Intel можно использовать инструкцию ассемблера fsincos или инструкцию expi (доступную, например, в D) для вычисления сложных

${\ displaystyle \ operatorname {expi} (z) = e ^ {iz} = \ cos (z) + i \ sin (z), \,}$

и просто разделите реальную и мнимую части.

Примечание. Чтобы явно вычислить комплексно-полярную форму, используйте следующие замены в общей форме:

Пусть и тогда ${\ Displaystyle г = {\ sqrt {-2 \ ln (и_ {1})}}}$${\ displaystyle z = 2 \ pi u_ {2}.}$

${\ displaystyle \ re ^ {iz} = {\ sqrt {-2 \ ln (u_ {1})}} e ^ {i2 \ pi u_ {2}} = {\ sqrt {-2 \ ln (u_ {1 })}} \ left [\ cos (2 \ pi u_ {2}) + i \ sin (2 \ pi u_ {2}) \ right].}$

Напротив, полярный метод здесь устраняет необходимость в вычислении косинуса и синуса. Вместо этого, решая точку на единичном круге, эти две функции можно заменить координатами x и y, нормированными на радиус. В частности, случайная точка ( x,  y) внутри единичной окружности проецируется на единичную окружность путем установки и формирования точки ${\ displaystyle {\ sqrt {x ^ {2} + y ^ {2}}}}$${\ displaystyle s = x ^ {2} + y ^ {2}}$

${\ displaystyle \ left ({\ frac {x} {\ sqrt {s}}}, {\ frac {y} {\ sqrt {s}}} \ right), \,}$

что является более быстрой процедурой, чем вычисление косинуса и синуса. Некоторые исследователи утверждают, что условная инструкция if (для отклонения точки за пределами единичного круга) может замедлить выполнение программ на современных процессорах, оснащенных конвейерной обработкой и предсказанием ветвлений. Также эта процедура требует примерно на 27% больше оценок основного генератора случайных чисел (только сгенерированные точки лежат внутри единичного круга). ${\ displaystyle \ pi / 4 \ приблизительно 79 \%}$

Затем эта случайная точка на окружности проецируется радиально на необходимое случайное расстояние с помощью

${\ Displaystyle {\ sqrt {-2 \ ln (s)}}, \,}$

с использованием тех же S, потому, что ей не зависит от случайной точки на окружности и сама равномерно распределяется от 0 до 1.

Простая реализация на Java с использованием среднего и стандартного отклонения:

private static double spare; private static boolean hasSpare = false; public static synchronized double generateGaussian(double mean, double stdDev) { if (hasSpare) { hasSpare = false; return spare * stdDev + mean; } else { double u, v, s; do { u = Math.random() * 2 - 1; v = Math.random() * 2 - 1; s = u * u + v * v; } while (s gt;= 1 || s == 0); s = Math.sqrt(-2.0 * Math.log(s) / s); spare = v * s; hasSpare = true; return mean + stdDev * u * s; } }

Не- поточно реализация в C ++, используя среднее и стандартное отклонение:

double generateGaussian(double mean, double stdDev) { static double spare; static bool hasSpare = false; if (hasSpare) { hasSpare = false; return spare * stdDev + mean; } else { double u, v, s; do { u = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0; v = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0; s = u * u + v * v; } while (s gt;= 1.0 || s == 0.0); s = sqrt(-2.0 * log(s) / s); spare = v * s; hasSpare = true; return mean + stdDev * u * s; } }

C ++ 11 GNU GCC libstdС ++ реализация std:: normal_distribution использует полярный метод Марсальи, как цитируется здесь.