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

Статистическое преобразование

Визуализация преобразования Бокса-Мюллера — цветные точки в единичном квадрате (u ₁ , u ₂ ), нарисованные в виде кругов, сопоставляются с двумерным гауссианом (z ₀ , z ₁ ), нарисованным крестиками. Графики на полях представляют собой функции распределения вероятностей z0 и z1. z0 и z1 неограничены; они кажутся в [-2,5, 2,5] из-за выбора иллюстрированных точек. В файле SVG наведите указатель мыши на точку, чтобы выделить ее и соответствующую ей точку.

Преобразование Бокса -Мюллера , разработанное Джорджем Эдвардом Пелэмом Боксом и Мервином Эдгаром Мюллером ^[1] , представляет собой метод выборки случайных чисел для генерации пар независимых , стандартных, нормально распределенных (нулевое ожидание , единичная дисперсия ) случайных чисел, учитывая источник равномерно распределенные случайные числа. Фактически, этот метод впервые был упомянут Раймондом Э.А.К. Пейли и Норбертом Винером в 1934 году ^{. [2]}

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

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

Основная форма

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

$${\displaystyle Z_{0}=R\cos(\Theta )={\sqrt {-2\ln U_{1}}}\cos(2\pi U_{2})\,}$$

$${\displaystyle Z_{1}=R\sin(\Theta )={\sqrt {-2\ln U_{1}}}\sin(2\pi U_{2}).\,}$$

Тогда Z0 и Z1 — _{независимые случайные} величины со стандартным нормальным распределением .

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

$${\displaystyle R^{2}=-2\cdot \ln U_{1}\,}$$

$${\displaystyle \Theta =2\pi U_{2}.\,}$$

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

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

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

Полярная форма была впервые предложена Дж. Беллом ^[6] и затем модифицирована Р. Кнопом. ^[7] Хотя было описано несколько различных версий полярного метода, здесь будет описана версия Р. Кнопа, поскольку она наиболее широко используется, отчасти из-за ее включения в числовые рецепты . Немного другая форма описана Д. Кнутом как «Алгоритм P» в книге «Искусство компьютерного программирования» . ^[8]

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

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

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

$${\displaystyle 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}}}}$$

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

Противопоставление двух форм

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

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

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

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

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

$${\displaystyle re^{iz}={\sqrt {-\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 деления для каждой нормальной переменной. В результате одно умножение и одна тригонометрическая функция заменяются одним делением и условным циклом.

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

Когда компьютер используется для создания однородной случайной величины, он неизбежно будет иметь некоторые неточности, поскольку существует нижняя граница того, насколько числа могут быть близки к 0. Если генератор использует 32 бита на выходное значение, это наименьшее ненулевое число, которое может быть создан . Когда и равны этому, преобразование Бокса-Мюллера дает нормальное случайное отклонение, равное . Это означает, что алгоритм не будет генерировать случайные величины со стандартными отклонениями более 6,660 от среднего значения. Это соответствует доле потерь из-за усечения, где – стандартное кумулятивное нормальное распределение. При использовании 64-битной версии предел смещается до стандартных отклонений, для которых . ${\displaystyle 2^{-32}}$ ${\displaystyle U_{1}}$ ${\displaystyle U_{2}}$ ${\textstyle \delta ={\sqrt {-2\ln(2^{-32})}}\cos(2\pi 2^{-32})\approx 6.660}$ ${\displaystyle 2(1-\Phi (\delta ))\simeq 2.738\times 10^{-11}}$ ${\displaystyle \Phi (\delta )}$ ${\displaystyle \delta =9.419}$ ${\displaystyle 2(1-\Phi (\delta ))<5\times 10^{-21}}$

Выполнение

С++

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

#include <cmath> #include <limits> #include <random> #include <utility>    // «мю» — среднее значение распределения, а «сигма» — стандартное отклонение. std :: pair < double , double > generateGaussianNoise ( double mu , double sigma ) { constexpr double two_pi = 2.0 * M_PI ;             //инициализируем генератор случайных унифицированных чисел (runif) в диапазоне от 0 до 1 static std :: mt19937 rng ( std :: random_device {}()); // Стандартный mersenne_twister_engine, заполненный с помощью rd() static std :: uniform_real_distribution <> runif ( 0.0 , 1.0 );         //создаем два случайных числа, убеждаемся, что u1 больше нуля double u1 , u2 ; do { u1 = runif ( rng ); } Пока ( u1 == 0 ); u2 = runif ( rng );                 //вычисление z0 и z1 auto mag = sigma * sqrt ( -2.0 * log ( u1 )); auto z0 = mag * cos ( two_pi * u2 ) + mu ; auto z1 = mag * sin ( two_pi * u2 ) + mu ;                             return std :: make_pair ( z0 , z1 ); }  

Юлия

"""  образец коробкимюллера(N)Сгенерируйте выборки `2N` из стандартного нормального распределения с помощью метода Бокса-Мюллера. """ функция boxmullersample ( N ) z = Array { Float64 }( undef , N , 2 ); for i в осях ( z , 1 ) z [ i , : ] .= sincospi ( 2 * rand ()); z [ i , : ] .*= sqrt ( - 2 * log ( rand ())); end vec ( z ) end                    """  boxmullersample(n,μ,σ)Сгенерируйте выборки n из нормального распределения со средним значением μ и стандартным отклонением σ, используя метод Бокса-Мюллера. """ функция boxmullersample ( n , μ , σ ) μ .+ σ * boxmullersample ( cld ( n , 2 ))[ 1 : n ]; end    

Смотрите также

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

Рекомендации

1. Коробка, GEP; Мюллер, Мервин Э. (1958). «Заметки о генерации случайных нормальных отклонений». Анналы математической статистики . 29 (2): 610–611. дои : 10.1214/aoms/1177706645 . JSTOR  2237361.
2. Раймонд Э.А.К. Пейли и Норберт Винер Преобразования Фурье в комплексной области, Нью-Йорк: Американское математическое общество (1934) §37.
3. Клоден и Платен, Численные решения стохастических дифференциальных уравнений , стр. 11–12.
4. Хоуз, Ли; Томас, Дэвид (2008). GPU Gems 3 — Эффективная генерация случайных чисел и применение с использованием CUDA . Pearson Education, Inc. ISBN 978-0-321-51526-1.
5. Шелдон Росс, Первый курс теории вероятностей , (2002), стр. 279–281.
6. Белл, Джеймс Р. (1968). «Алгоритм 334: Нормальные случайные отклонения». Коммуникации АКМ . 11 (7): 498. дои : 10.1145/363397.363547 .
7. Кноп, Р. (1969). «Замечание к алгоритму 334 [G5]: Нормальные случайные отклонения». Коммуникации АКМ . 12 (5): 281. дои : 10.1145/362946.362996 .
8. Кнут, Дональд (1998). Искусство компьютерного программирования: Том 2: Получисловые алгоритмы . п. 122. ИСБН 0-201-89684-2.
9. Эверетт Ф. Картер-младший, Генерация и применение случайных чисел, Четвертые измерения (1994), Vol. 16, № 1 и 2.
10. Оценка 2 π U ₁ считается за одно умножение, поскольку значение 2 π можно вычислить заранее и использовать повторно.

Внешние ссылки

• Вайсштейн, Эрик В. «Трансформация Бокса-Мюллера». Математический мир .
• Как преобразовать равномерное распределение в распределение Гаусса (код C)