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

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

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

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

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

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

Предположим, что U ₁ и U ₂ являются независимыми выборками, выбранными из равномерного распределения на единичном интервале (0, 1) . Позволять $$Z_{0}=R\cos(\Theta )={\sqrt {-2\ln U_{1}}}\cos(2\pi U_{2})\,$$ и $$Z_{1}=R\sin(\Theta )={\sqrt {-2\ln U_{1}}}\sin(2\pi U_{2}).\,$$

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

Вывод основан на свойстве двумерной декартовой системы , где координаты X и Y описываются двумя независимыми и нормально распределенными случайными величинами, случайными величинами для R ² и Θ (показано выше) в соответствующих полярных координатах также независимы и могут быть выражены как $$R^{2}=-2\cdot \ln U_{1}\,$$ и $$\Theta =2\pi U_{2}.\,$$

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

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

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

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

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

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

Теперь мы отождествляем значение s со значением U ₁ и $\theta /(2\pi )$ с U ₂ в основной форме. Как показано на рисунке, значения $\cos \theta =\cos 2\pi U_{2}$ и $\sin \theta =\sin 2\pi U_{2}$ в основной форме можно заменить соотношениями $\cos \theta =u/R=u/{\sqrt {s}}$ и ${\textstyle \sin \theta =v/R=v/{\sqrt {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}={\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)=e^{iz}=\cos z+i\sin z,\,$$ и просто разделите действительную и мнимую части.

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

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

Выполнение

С++

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

#include <cmath>
#include <limits>
#include <random>
#include <utility>

//"mu" is the mean of the distribution, and "sigma" is the standard deviation.
std::pair<double, double> generateGaussianNoise(double mu, double sigma)
{
    constexpr double two_pi = 2.0 * M_PI;

    //initialize the random uniform number generator (runif) in a range 0 to 1
    static std::mt19937 rng(std::random_device{}()); // Standard mersenne_twister_engine seeded with rd()
    static std::uniform_real_distribution<> runif(0.0, 1.0);

    //create two random numbers, make sure u1 is greater than zero
    double u1, u2;
    do
    {
        u1 = runif(rng);
    }
    while (u1 == 0);
    u2 = runif(rng);

    //compute z0 and 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);
}

JavaScript

function rand_normal()  {

    /* Syntax:
     *
     *    [ x, y ]  =   rand_normal();
     *      x       =   rand_normal()[0];
     *      y       =   rand_normal()[1];
     *
    */


    // Box-Muller Transform:

    let phi  =   2 * Math.PI * Math.random();
    let R    =   Math.sqrt(  -2 * Math.log( Math.random() )  );
    let x    =   R * Math.cos(phi);
    let y    =   R * Math.sin(phi);


    // Return values:

    return [ x, y ];

}

Юлия

"""
    boxmullersample(N)

Generate `2N` samples from the standard normal distribution using the Box-Muller method.
"""
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,μ,σ)

Generate `n` samples from the normal distribution with mean `μ` and standard deviation `σ` using the Box-Muller method.
"""
function boxmullersample(n,μ,σ)
    μ .+ σ*boxmullersample(cld(n,2))[1:n];
end

См. также

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

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

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