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

Марсальи . Полярный метод — это метод выборки псевдослучайных чисел для генерации пары независимых стандартных нормальных случайных величин .

Стандартные нормальные случайные величины часто используются в информатике , вычислительной статистике и, в частности, в приложениях метода Монте-Карло .

Полярный метод работает путем выбора случайных точек ( x , y ) в квадрате −1 < x < 1, −1 < y < 1 до тех пор, пока

$0<s=x^{2}+y^{2}<1,\,$

а затем возвращая требуемую пару нормальных случайных величин как

$x{\sqrt {\frac {-2\ln(s)}{s}}}\,,\ \ y{\sqrt {\frac {-2\ln(s)}{s}}},$

или, что то же самое,

${\frac {x}{\sqrt {s}}}{\sqrt {-2\ln(s)}}\,,\ \ {\frac {y}{\sqrt {s}}}{\sqrt {-2\ln(s)}},$

где $x/{\sqrt {s}}$ и $y/{\sqrt {s}}$ представляют косинус и синус угла, который вектор ( x , y ) образует с x осью .

Теоретическая основа

Основополагающую теорию можно резюмировать следующим образом:

Если u равномерно распределено в интервале 0 ≤ u < 1, то точка (cos(2π u ), sin(2π u )) равномерно распределен по единичной окружности х ² + и ² = 1, и умножив эту точку на независимую случайная величина ρ, распределение которой равно

$\Pr(\rho <a)=\int _{0}^{a}re^{-r^{2}/2}\,dr$

поставит точку

$\left(\rho \cos(2\pi u),\rho \sin(2\pi u)\right)$

координаты которых совместно распределены как два независимых стандарта обычные случайные величины.

История

Эта идея восходит к Лапласу , которому Гаусс приписывает открытие вышеизложенного.

$I=\int _{-\infty }^{\infty }e^{-x^{2}/2}\,dx$

извлекая квадратный корень из

$I^{2}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }e^{-(x^{2}+y^{2})/2}\,dx\,dy=\int _{0}^{2\pi }\int _{0}^{\infty }re^{-r^{2}/2}\,dr\,d\theta .$

Преобразование в полярные координаты показывает, что θ есть равномерно распределены (постоянная плотность) от 0 до 2π, и что радиальное расстояние r имеет плотность

$re^{-r^{2}/2}.\,$

( р ² имеет соответствующее распределение хи-квадрат .)

Этот метод создания пары независимых стандартных нормальных переменных путем радиального проецирования случайной точки на единичную окружность на расстояние, определяемое квадратным корнем из переменной хи-квадрат-2, называется полярным методом генерации пары нормальных случайных величин. ,

Практические соображения

Прямое применение этой идеи

$x={\sqrt {-2\ln(u_{1})}}\cos(2\pi u_{2}),\quad y={\sqrt {-2\ln(u_{1})}}\sin(2\pi u_{2})$

называется преобразованием Бокса-Мюллера , в котором переменная хи обычно равна генерируется как

${\sqrt {-2\ln(u_{1})}},$

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

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

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

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

Позволять $r={\sqrt {-2\ln(u_{1})}}$ и $z=2\pi u_{2}.$ Затем

$\ 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].$

Напротив, полярный метод здесь избавляет от необходимости вычислять косинус и синус. Вместо этого, найдя точку на единичном круге, эти две функции можно заменить координатами x и y, нормализованными к ${\sqrt {x^{2}+y^{2}}}$ радиус. В частности, случайная точка ( x , y ) внутри единичного круга проецируется на единичную окружность, установив $s=x^{2}+y^{2}$ и формируем точку

$\left({\frac {x}{\sqrt {s}}},{\frac {y}{\sqrt {s}}}\right),\,$

это более быстрая процедура, чем вычисление косинуса и синуса. Некоторые исследователи утверждают, что условная инструкция if (для отклонения точки за пределами единичного круга) может замедлить работу программ на современных процессорах, оснащенных конвейерной обработкой и прогнозированием ветвей. Кроме того, эта процедура требует примерно на 27% больше вычислений базового генератора случайных чисел (только $\pi /4\approx 79\%$ сгенерированных точек лежат внутри единичного круга).

Затем эта случайная точка на окружности проецируется радиально на необходимое случайное расстояние с помощью

${\sqrt {-2\ln(s)}},\,$

используя одно и то же s , поскольку оно не зависит от случайной точки на окружности и само по себе равномерно распределено от 0 до 1.

Выполнение

Ява

Простая реализация на Java с использованием среднего и стандартного отклонения:

private static double spare;
private static boolean hasSpare = false;

public static synchronized double generateGaussian(double mean, double stdDev) {
    if (hasSpare) {
        hasSpare = false;
        return spare * stdDev + mean;
    } else {
        double u, v, s;
        do {
            u = Math.random() * 2 - 1;
            v = Math.random() * 2 - 1;
            s = u * u + v * v;
        } while (s >= 1 || s == 0);
        s = Math.sqrt(-2.0 * Math.log(s) / s);
        spare = v * s;
        hasSpare = true;
        return mean + stdDev * u * s;
    }
}

С++

Непотокобезопасная с реализация на C++ использованием среднего и стандартного отклонения:

double generateGaussian(double mean, double stdDev) {
    static double spare;
    static bool hasSpare = false;

    if (hasSpare) {
        hasSpare = false;
        return spare * stdDev + mean;
    } else {
        double u, v, s;
        do {
            u = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;
            v = (rand() / ((double)RAND_MAX)) * 2.0 - 1.0;
            s = u * u + v * v;
        } while (s >= 1.0 || s == 0.0);
        s = sqrt(-2.0 * log(s) / s);
        spare = v * s;
        hasSpare = true;
        return mean + stdDev * u * s;
    }
}

в C++11 в GNU GCC libstdc++ Реализация std::normal_distribution использует полярный метод Марсальи, как указано здесь .

Юлия

Простая реализация Джулии :

"""
    marsagliasample(N)

Generate `2N` samples from the standard normal distribution using the Marsaglia method.
"""
function marsagliasample(N)
    z = Array{Float64}(undef,N,2);
    for i in axes(z,1)
        s = Inf;
        while s > 1
            z[i,:] .= 2rand(2) .- 1;
            s = sum(abs2.(z[i,:]))
        end
        z[i,:] .*= sqrt(-2log(s)/s);
    end
    vec(z)
end

"""
    marsagliasample(n,μ,σ)

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

Цикл for можно распараллелить с помощью Threads.@threads макрос.