Коробка -маллер преобразование

Визуализация преобразования коробки -маллер - цветные точки на квадрате единицы (U ₁ , U ₂ ), нарисованные в виде кругов, отображаются с 2 -й гауссовой (Z ₀ , Z ₁ ), нарисованной в виде скрещиваний. Графики на полях - это функции распределения вероятностей Z0 и Z1. Z0 и Z1 не ограничены; По -видимому, они находятся в [-2,5, 2,5] из -за выбора иллюстрированных точек. В файле SVG парят над точкой, чтобы выделить его и его соответствующую точку.

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

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

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

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

Предположим, что 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}).\,$$

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

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

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

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

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

Два равномерно распределенных значения, U и V используются для получения значения S = R ² , который также равномерно распределен. Определения синуса и косинуса затем применяются к основной форме преобразования коробки -маллинг, чтобы избежать использования тригонометрических функций.

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

Учитывая U и V , независимое и равномерно распределено в замкнутом интервале [-1, +1] , установить s = r ² = u ² + 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% от общего входа, сгенерированных равномерно распределенных распределенных чисел, IE отбрасывает 4/ π - 1 ≈ 27,32% равномерно распределенные пары случайных чисел на гауссовые пары случайных чисел, сгенерированная, требуя 4/ π ≈ 1,2732 Случайные числа на вывод случайное число.

Основная форма требует двух умножений, 1/2 логарифма, 1/2 квадратного корня и одной тригонометрической функции для каждого нормального вариата. На некоторых процессорах косинус и синус того же аргумента могут быть рассчитаны параллельно с помощью одной инструкции. Примечательно, что для машин на основе Intel можно использовать инструкцию Fsincos Assembler или инструкцию 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}$.

Выполнение

C ++

Стандартное преобразование коробки - маллинг генерирует значения из стандартного нормального распределения ( т.е. стандартные нормальные отклонения ) со средним 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

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

• Обратная выборка преобразования
• Marsaglia Polar Method , аналогичный преобразованию с Box -Muller, который использует картезианские координаты вместо полярных координат

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

• Вейсштейн, Эрик У. "Трансформация боксов" . MathWorld .
• Как преобразовать равномерное распределение в гауссовое распределение (C код)