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

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

Визуализация преобразования Бокса–Мюллера — цветные точки в единичном квадрате (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}).\,}$$

Тогда Z ₀ и Z ₁ — независимые случайные величины со стандартным нормальным распределением .

Вывод ^[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] Хотя было описано несколько различных версий полярного метода, здесь будет описана версия Р. Кнопа, поскольку она наиболее широко используется, отчасти из-за ее включения в Numerical Recipes . Немного иная форма описана как «Алгоритм P» Д. Кнутом в The Art of Computer Programming . ^[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>    //"mu" - это среднее значение распределения, а "sigma" - это стандартное отклонение. 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 ); } while ( 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 ;                             возврат std :: make_pair ( z0 , z1 ); }  

JavaScript

функция rand_normal () {   /* Синтаксис:  *  * [ x, y ] = rand_normal();  * x = rand_normal()[0];  * y = rand_normal()[1];  *  */ // Преобразование Бокса-Мюллера: пусть phi = 2 * Math . PI * Math . random (); пусть R = Math . sqrt ( - 2 * Math . log ( Math . random () ) ); пусть x = R * Math . cos ( phi ); пусть y = R * Math . sin ( phi );                              // Возвращаемые значения: вернуть [ х , у ];    }

Джулия

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

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

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

Ссылки

1. Бокс, GEP; Мюллер, Мервин Э. (1958). «Заметка о генерации случайных нормальных отклонений». Анналы математической статистики . 29 (2): 610–611. doi : 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. Bell, James R. (1968). "Алгоритм 334: Нормальные случайные отклонения". Сообщения ACM . 11 (7): 498. doi : 10.1145/363397.363547 .
7. Кноп, Р. (1969). «Замечание об алгоритме 334 [G5]: Нормальные случайные отклонения». Сообщения ACM . 12 (5): 281. doi : 10.1145/362946.362996 .
8. Кнут, Дональд (1998). Искусство программирования: Том 2: Получисленные алгоритмы . стр. 122. ISBN 0-201-89684-2.
9. Эверетт Ф. Картер-младший, Генерация и применение случайных чисел, Forth Dimensions (1994), том 16, № 1 и 2.
10. Оценка 2 π U ₁ считается одним умножением, поскольку значение 2 π можно вычислить заранее и использовать повторно.

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

• Вайсштейн, Эрик В. «Преобразование Бокса-Мюллера». MathWorld .
• Как преобразовать равномерное распределение в гауссово (код C)