Преобразование Бокса – Мюллера

редактировать

Визуализация преобразования Бокса – Мюллера - цветные точки в единичном квадрате (u1, u2), нарисованные в виде кругов, отображаются в 2D-гауссиане (z0, z1), нарисованы крестами. Графики на полях - это функции распределения вероятностей z0 и z1. Отметим, что z0 и z1 неограничены; они появляются в [-2,5,2,5] из-за выбора проиллюстрированных точек. В файле SVG наведите указатель мыши на точку, чтобы выделить ее и соответствующую ей точку.

Преобразование Бокса – Мюллера, автор Джордж Эдвард Пелхэм Бокс и Мервин Эдгар Мюллер, это метод выборки случайных чисел для генерации пар независимых, стандартных, нормально распределенных (нулевое ожидание, единица дисперсия ) случайных чисел, заданных источником равномерно распределенных случайных чисел. Фактически, этот метод впервые был упомянут в явной форме Раймондом Э. А. С. Пэли и Норбертом Винером в 1934 году.

Преобразование Бокса – Мюллера обычно выражается в двух формах. Основная форма, данная Боксом и Мюллером, берет две выборки из равномерного распределения на интервале [0, 1] и отображает их на две стандартные, нормально распределенные выборки. Полярная форма берет две выборки из другого интервала, [-1, +1], и отображает их в две нормально распределенные выборки без использования синусоидальных или косинусных функций.

Преобразование Бокса – Мюллера было разработано как более эффективная с вычислительной точки зрения альтернатива методу выборки обратного преобразования. Алгоритм зиккурата дает более эффективный метод для скалярных процессоров (например, старых ЦП), в то время как преобразование Бокса – Мюллера лучше для процессоров с векторными модулями (например, графические процессоры или современные процессоры).

Содержание

1 Базовая форма
2 Полярная форма
3 Противопоставление двух форм
4 Усечение хвостов
5 Реализация
6 См. Также
7 Ссылки
8 Внешние ссылки

Базовые form

Предположим, что U 1 и U 2 являются независимыми выборками, выбранными из равномерного распределения на единичном интервале (0, 1). Пусть

Z 0 знак равно R соз ⁡ (Θ) = - 2 пер ⁡ U 1 соз ⁡ (2 π U 2) {\ displaystyle Z_ {0} = R \ cos (\ Theta) = {\ sqrt {-2 \ ln U_ {1}}} \ cos (2 \ pi U_ {2}) \,}

Z _0 = R \ соз (\ Theta) = \ sqrt {-2 \ ln U_1} \ cos (2 \ pi U_2) \,

Z 1 = R sin ⁡ (Θ) = - 2 ln ⁡ U 1 sin ⁡ (2 π U 2). {\ displaystyle Z_ {1} = R \ sin (\ Theta) = {\ sqrt {-2 \ ln U_ {1}}} \ sin (2 \ pi U_ {2}). \,}

Z_1 = R \ sin (\ Theta) = \ sqrt {-2 \ ln U_1} \ sin (2 \ pi U_2). \,

Тогда Z 0 и Z 1 являются независимыми случайными величинами со стандартным нормальным распределением .

Вывод основан на свойстве двумерного Декартова система, где координаты X и Y описываются двумя независимыми и нормально распределенными случайными величинами, случайные величины для R и Θ (показанные выше) в соответствующих полярных координатах также независимы и могут быть выражены как

R 2 = - 2 ⋅ пер U 1 {\ displaystyle R ^ {2} = - 2 \ cdot \ ln U_ {1} \,}

R ^ 2 = -2 \ cdot \ ln U_1 \,

Θ = 2 π U 2. {\ displaystyle \ Theta = 2 \ pi U_ {2}. \,}

\ Theta = 2 \ pi U_2. \,

Поскольку R является квадратом нормы стандартной двумерной нормальной переменной (X, Y), она имеет распределение хи-квадрат с двумя степенями свободы. В частном случае двух степеней свободы распределение хи-квадрат совпадает с экспоненциальным распределением , а приведенное выше уравнение для R представляет собой простой способ генерировать требуемую экспоненциальную переменную.

Полярная форма

Два равномерно распределенных значения u и v используются для получения значения s = R, которое также равномерно распределено. Затем определения синуса и косинуса применяются к базовой форме преобразования Бокса – Мюллера, чтобы избежать использования тригонометрических функций.

Полярная форма была сначала предложена Дж. Беллом, а затем модифицирована Р. Кнопом. Хотя было описано несколько различных версий полярного метода, здесь будет описана версия Р. Кнопа, потому что она наиболее широко используется, отчасти из-за ее включения в Числовые рецепты.

Учитывая u и v, независимые и равномерно распределенные в закрытом интервале [−1, +1], положим s = R = u + v. Если s = 0 или s ≥ 1, отбросьте u и v и попробуйте другую пару (u, v). Поскольку u и v равномерно распределены и допускаются только точки внутри единичного круга, значения s будут равномерно распределены и в открытом интервале (0, 1). Последнее можно увидеть, вычислив кумулятивную функцию распределения для s в интервале (0, 1). Это площадь круга радиусом $s {\ displaystyle \ scriptstyle {\ sqrt {s}}}$ $\ scriptstyle \ sqrt {s}$ , деленная на $π {\ displaystyle \ scriptstyle \ pi}$ $\ scriptstyle \ пи$ . Отсюда мы находим, что функция плотности вероятности имеет постоянное значение 1 на интервале (0, 1). Точно так же угол θ, деленный на $2 π {\ displaystyle \ scriptstyle 2 \ pi}$ $\ scriptstyle 2 \ pi$ , равномерно распределен в интервале [0, 1) и не зависит от s.

Теперь мы идентифицируем значение s со значением U 1 и $θ / (2 π) {\ displaystyle \ scriptstyle \ theta / (2 \ pi)}$ $\ scriptstyle \ theta / (2 \ pi)$ с U 2 в основной форме. Как показано на рисунке, значения $cos ⁡ θ = cos ⁡ 2 π U 2 {\ displaystyle \ scriptstyle \ cos \ theta = \ cos 2 \ pi U_ {2}}$ $\ scriptstyle \ cos \ theta = \ cos 2 \ pi U_2$ и $грех ⁡ θ = грех ⁡ 2 π U 2 {\ displaystyle \ scriptstyle \ sin \ theta = \ sin 2 \ pi U_ {2}}$ $\ scriptstyle \ sin \ theta = \ sin 2 \ pi U_2$ в базовой форме может быть заменен на коэффициенты $соз ⁡ θ = u / R = u / s {\ displaystyle \ scriptstyle \ cos \ theta = u / R = u / {\ sqrt {s}}}$ $\ scriptstyle \ cos \ theta = u / R = u / \ sqrt {s}$ и $грех ⁡ θ = v / R = v / s {\ displaystyle \ scriptstyle \ sin \ theta = v / R = v / {\ sqrt {s}}}$ $\ scriptstyle \ sin \ theta = v / R = v / \ sqrt {s}$ соответственно. Преимущество состоит в том, что можно избежать прямого вычисления тригонометрических функций. Это полезно, когда вычисление тригонометрических функций обходится дороже, чем одно деление, которое заменяет каждую из них.

Так же, как базовая форма дает два стандартных нормальных отклонения, так и этот альтернативный расчет.

z 0 знак равно - 2 пер U 1 соз ⁡ (2 π U 2) = - 2 пер ⁡ s (нас) = U ⋅ - 2 пер ⁡ сс {\ displaystyle z_ {0} = {\ sqrt {- 2 \ ln U_ {1}}} \ cos (2 \ pi U_ {2}) = {\ sqrt {-2 \ ln s}} \ left ({\ frac {u} {\ sqrt {s}}} \ справа) = u \ cdot {\ sqrt {\ frac {-2 \ ln s} {s}}}}

z_0 = \ sqrt {-2 \ ln U_1} \ cos (2 \ pi U_2) = \ sqrt {-2 \ ln s} \ left (\ frac {u} {\ sqrt {s}} \ right) = u \ cdot \ sqrt {\ frac {-2 \ ln s} {s}}

z 1 = - 2 ln ⁡ U 1 sin ⁡ (2 π U 2) = - 2 ln ⁡ s (vs) = v ⋅ - 2 ln ⁡ ss. {\ displaystyle z_ {1} = {\ sqrt {-2 \ ln U_ {1}}} \ sin (2 \ pi U_ {2}) = {\ sqrt {-2 \ ln s}} \ left ({\ frac {v} {\ sqrt {s}}} \ right) = v \ cdot {\ sqrt {\ frac {-2 \ ln s} {s}}}.}

z_1 = \ sqrt {-2 \ ln U_1} \ sin (2 \ pi U_2) = \ sqrt {-2 \ ln s} \ left (\ frac {v} {\ sqrt {s} } \ right) = v \ cdot \ sqrt {\ frac {-2 \ ln s} {s}}.

Сравнение двух форм

Полярный метод отличается от базового метода тем, что он представляет собой тип отбраковки. Он отбрасывает некоторые сгенерированные случайные числа, но может быть быстрее, чем базовый метод, потому что его проще вычислить (при условии, что генератор случайных чисел относительно быстр) и он более устойчив в числовом отношении. Избегание использования дорогих тригонометрических функций увеличивает скорость по сравнению с базовой формой. Он отбрасывает 1 - π / 4 ≈ 21,46% от общего числа сгенерированных пар случайных чисел с равномерным распределением входных данных, то есть отбрасывает 4 / π - 1 ≈ 27,32% равномерно распределенных пар случайных чисел на каждую сгенерированную пару случайных чисел гауссова, что требует 4 / π ≈ 1,2732 входных случайных числа на каждое выходное случайное число.

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

exp ⁡ (iz) = eiz = cos ⁡ z + i sin ⁡ z, {\ displaystyle \ exp (iz) = e ^ {iz} = \ cos z + i \ sin z, \,}

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

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

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

Пусть $r = - 2 ln ⁡ (u 1) {\ displaystyle r = {\ sqrt {-2 \ ln (u_ {1})}}}$ ${\ displaystyle r = {\ sqrt {-2 \ ln (u_ {1})}}}$ и $z = 2 π u 2. {\ displaystyle z = 2 \ pi u_ {2}.}$ ${\ displaystyle z = 2 \ pi u_ {2}.}$ Тогда

reiz = - 2 ln ⁡ (u 1) ei 2 π u 2 = - 2 ln ⁡ (u 1) [cos ⁡ (2 π u 2) + я sin ⁡ (2 π 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].}

{\ 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].}

Полярная форма требует умножения на 3/2, 1/2 логарифм, 1/2 квадратного корня и 1/2 деления для каждой нормальной вариации. Эффект заключается в замене одного умножения и одной тригонометрической функции на одно деление и условный цикл.

Усечение хвостов

Когда компьютер используется для создания однородной случайной величины, он неизбежно будет иметь некоторые неточности, потому что существует нижняя граница того, насколько близки числа к нулю. Если генератор использует 32 бита на выходное значение, наименьшее ненулевое число, которое может быть сгенерировано: $2–32 {\ displaystyle 2 ^ {- 32}}$ $2 ^ {- 32}$ . Когда $U 1 {\ displaystyle U_ {1}}$ $U_ {1}$ и $U 2 {\ displaystyle U_ {2}}$ $U_ {2}$ равны этому, преобразование Бокса – Мюллера дает нормальное случайное отклонение, равное $δ = - 2 ln ⁡ (2–32) cos ⁡ (2 π 2–32) ≈ 6.660 {\ displaystyle \ delta = {\ sqrt {-2 \ ln (2 ^ {- 32})}} \ cos (2 \ pi 2 ^ {- 32}) \ приблизительно 6,660}$ ${\ displaystyle \ delta = {\ sqrt {-2 \ ln (2 ^ {- 32})}} \ cos (2 \ pi 2 ^ {- 32 }) \ приблизительно 6,660}$ . Это означает, что алгоритм не будет генерировать случайные величины, превышающие 6,660 стандартных отклонений от среднего. Это соответствует пропорции $2 (1 - Φ (δ)) ≃ 2.738 × 10 - 11 {\ displaystyle 2 (1- \ Phi (\ delta)) \ simeq 2.738 \ times 10 ^ {- 11}}$ ${\ displaystyle 2 (1- \ Phi (\ delta)) \ simeq 2.738 \ times 10 ^ {- 11}}$ потеряны из-за усечения, где $Φ (δ) {\ displaystyle \ Phi (\ delta)}$ ${\ displaystyle \ Phi (\ delta)}$ - стандартное кумулятивное нормальное распределение. С 64 битами предел сдвигается до $δ = 9,419 {\ displaystyle \ delta = 9,419}$ ${\ displaystyle \ delta = 9,419}$ стандартных отклонений, для которых $2 (1 - Φ (δ)) < 5 × 10 − 21 {\displaystyle 2(1-\Phi (\delta))<5\times 10^{-21}}$ ${\ displaystyle 2 (1- \ Phi (\ delta)) <5 \ раз 10 ^ {- 21}}$ .

Реализация

Стандартное преобразование Бокса – Мюллера генерирует значения из стандартного нормального распределения (т. Е. стандартных нормальных отклонений ) со средним 0 и стандартным отклонением 1. Реализация ниже в стандарте C ++ генерирует значения из любого нормального распределения со средним значением $μ {\ displaystyle \ mu}$ $\ mu$ и дисперсией $σ 2 {\ displaystyle \ sigma ^ {2}}$ $\ sigma ^ {2}$ . Если $Z {\ displaystyle Z}$ $Z$ является стандартным нормальным отклонением, то $X = Z σ + μ {\ displaystyle X = Z \ sigma + \ mu}$ $X = Z \ sigma + \ mu$ будет иметь нормальное распределение со средним значением $μ {\ displaystyle \ mu}$ $\ mu$ и стандартным отклонением $σ {\ displaystyle \ sigma}$ $\ sigma$ . Обратите внимание, что генератор случайных чисел был засеян, чтобы гарантировать, что новые псевдослучайные значения будут возвращены из последовательных вызовов функции generateGaussianNoise.

#include #include #include double generateGaussianNoise (двойной mu, двойной сигма) {constexpr double epsilon = std :: numeric_limits :: epsilon (); constexpr double two_pi = 2.0 * M_PI; статический std :: mt19937 rng (std :: random_device {} ()); // Стандартный mersenne_twister_engine, засеянный rd () static std :: uniform_real_distribution <>runif (0.0, 1.0); двойной u1, u2; сделать {u1 = runif (rng); u2 = runif (rng); } while (u1 <= epsilon); auto z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2); // auto z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2); // not used here! return z0 * sigma + mu; }

См. также

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

Ссылки

Howes, Ли; Томас, Дэвид (2008). GPU Gems 3 - Эффективное создание случайных чисел и применение с использованием CUDA. Pearson Education, Inc. ISBN 978-0-321-51526-1. CS1 maint: ref = harv (ссылка )

^Box, GEP; Muller, Mervin E. (1958). «Заметка о генерации случайных нормальных отклонений». Анналы математической статистики. 29 (2): 610–611. doi : 10.1214 / aoms / 1177706645. JSTOR 2237361.
^Реймонд Е.А. Пэли и Норберт Винер преобразования Фурье в комплексной области, Нью-Йорк: Американское математическое общество (1934) §37.
^Клоеден и Платен, Численные решения стохастических дифференциальных уравнений, стр. 11–12
^Howes Thomas 2008.
^Шелдон Росс, Первый курс теории вероятностей, (2002), стр. 279–281
^ B элл, Джеймс Р. (1968). «Алгоритм 334: нормальные случайные отклонения». Коммуникации ACM. 11 (7): 498. doi : 10.1145 / 363397.363547.
^Кноп, Р. (1969). «Замечание по алгоритму 334 [G5]: Нормальные случайные отклонения». Коммуникации ACM. 12 (5): 281. doi : 10.1145 / 362946.362996.
^Эверетт Ф. Картер младший, Генерация и применение случайных чисел, Forth Dimensions (1994), Т. 16, № 1 и 2.
^Обратите внимание, что оценка 2πU 1 засчитывается как одно умножение, потому что значение 2π может быть вычислено заранее и использоваться повторно.