Преобразование равномерно распределенной случайной величины в нормально распределенную

Этот вопрос уже давно подробно изучен, и наиболее широкое распространение получил метод полярных координат, предложенный Джорджем Боксом, Мервином Мюллером и Джорджем Марсальей в 1958 году. Данный метод позволяет получить пару независимых нормально распределенных случайных величин с математическим ожиданием 0 и дисперсией 1 следующим образом:



где Z₀ и Z₁ — искомые значения, s = u² + v², а u и v — равномерно распределенные на отрезке (-1, 1) случайные величины, подобранные таким образом, чтобы выполнялось условие 0

Многие используют эти формулы, даже не задумываясь, а многие даже и не подозревают об их существовании, так как пользуются готовыми реализациями. Но есть люди, у которых возникают вопросы: «Откуда взялась эта формула? И почему получается сразу пара величин?». Далее я постараюсь дать наглядный ответ на эти вопросы.

Для начала напомню, что такое плотность вероятности, функция распределения случайной величины и обратная функция. Допустим, есть некая случайная величина, распределение которой задано функцией плотности f(x), имеющей следующий вид:

Это означает, что вероятность того, что значение данной случайной величины окажется в интервале (A, B), равняется площади затененной области. И как следствие, площадь всей закрашенной области должна равняться единице, так как в любом случае значение случайной величины попадет в область определения функции f.

Функция распределения случайной величины является интегралом от функции плотности. И в данном случае ее примерный вид будет такой:

Тут смысл в том, что значение случайной величины будет меньше чем A с вероятностью B. И как следствие, функция никогда не убывает, а ее значения лежат в отрезке [0, 1].

Обратная функция — это функция, которая возвращает аргумент исходной функции, если в нее передать значение исходной функции. Например, для функции x² обратной будет функция извлечения корня, для sin(x) это arcsin(x) и т.д.

Так как большинство генераторов псевдослучайных чисел на выходе дают только равномерное распределение, то часто возникает необходимость его преобразования в какое-либо другое. В данном случае в нормальное Гауссовское:

Основу всех методов преобразования равномерного распределения в любое другое составляет метод обратного преобразования. Работает он следующим образом. Находится функция, обратная функции необходимого распределения, и в качестве аргумента передается в нее равномерно распределенная на отрезке (0, 1) случайная величина. На выходе получаем величину с требуемым распределением. Для наглядности привожу следующую картинку.

Таким образом, равномерный отрезок как бы размазывается в соответствии с новым распределением, проецируясь на другую ось через обратную функцию. Но проблема в том, что интеграл от плотности Гауссовского распределения вычисляется непросто, поэтому вышеперечисленным ученым пришлось схитрить.

Существует распределение хи-квадрат (распределение Пирсона), которое представляет собой распределение суммы квадратов k независимых нормальных случайных величин. И в случае, когда k = 2, это распределение является экспоненциальным.

Это означает, что если у точки в прямоугольной системе координат будут случайные координаты X и Y, распределенные нормально, то после перевода этих координат в полярную систему (r, θ) квадрат радиуса (расстояния от начала координат до точки) будет распределен по экспоненциальному закону, так как квадрат радиуса — это сумма квадратов координат (по закону Пифагора). Плотность распределения таких точек на плоскости будет выглядеть следующим образом:

Так как она равноценна во всех направлениях, угол θ будет иметь равномерное распределение в диапазоне от 0 до 2π. Справедливо и обратное: если задать точку в полярной системе координат с помощью двух независимых случайных величин (угла, распределенного равномерно, и радиуса, распределенного экспоненциально), то прямоугольные координаты этой точки будут являться независимыми нормальными случайными величинами. А экспоненциальное распределение из равномерного получить уже гораздо проще, с помощью того же метода обратного преобразования. В этом и заключается суть полярного метода Бокса-Мюллера.

Теперь выведем формулы.

(1)

Для получения r и θ нужно сгенерировать две равномерно распределенные на отрезке (0, 1) случайные величины (назовем их u и v), распределение одной из которых (допустим v) необходимо преобразовать в экспоненциальное для получения радиуса. Функция экспоненциального распределения выглядит следующим образом:

Обратная к ней функция:

Так как равномерное распределение симметрично, то аналогично преобразование будет работать и с функцией

Из формулы распределения хи-квадрат следует, что λ = 0,5. Подставим в эту функцию λ, v и получим квадрат радиуса, а затем и сам радиус:

Угол получим, растянув единичный отрезок до 2π:

Теперь подставим r и θ в формулы (1) и получим:

(2)

Эти формулы уже готовы к использованию. X и Y будут независимы и распределены нормально с дисперсией 1 и математическим ожиданием 0. Чтобы получить распределение с другими характеристиками достаточно умножить результат функции на среднеквадратическое отклонение и прибавить математическое ожидание.

Но есть возможность избавиться от тригонометрических функций, задав угол не прямо, а косвенно через прямоугольные координаты случайной точки в круге. Тогда через эти координаты можно будет вычислить длину радиус-вектора, а потом найти косинус и синус, поделив на нее x и y соответственно. Как и почему это работает?

Выберем случайную точку из равномерно распределенных в круге единичного радиуса и обозначим квадрат длины радиус-вектора этой точки буквой s:

Выбор осуществляется заданием случайных прямоугольных координат x и y, равномерно распределенных в интервале (-1, 1), и отбрасыванием точек, которые не принадлежат кругу, а также центральной точки, в которой угол радиус-вектора не определен. То есть должно выполниться условие 0 < s < 1. Тогда, как и в случае с Гауссовским распределением на плоскости, угол θ будет распределен равномерно. Это очевидно — количество точек в каждом направлении одинаково, значит каждый угол равновероятен. Но есть и менее очевидный факт — s тоже будет иметь равномерное распределение. Хоть s и является функцией от x и y, от каждой из этих координат в отдельности s будет независима. Поэтому мы можем воспользоваться значением s для получения экспоненциального распределения, не генерируя третью случайную величину. Подставим теперь s в формулы (2) вместо v, а вместо тригонометрических функций — их расчет делением координаты на длину радиус-вектора, которая в данном случае является корнем из s:

Получаем формулы, как в начале статьи. Недостаток этого метода — отбрасывание точек, не вошедших в круг. То есть использование только 78,5% сгенерированных случайных величин. На старых компьютерах отсутствие тригонометрических функций всё равно давало большое преимущество. Сейчас, когда одна команда процессора за мгновение вычисляет одновременно синус и косинус, думаю, эти методы могут еще посоревноваться.

Лично у меня остается еще два вопроса:

• Почему значение s распределено равномерно?


• Почему сумма квадратов двух нормальных случайных величин распределена экспоненциально?

Так как s — это квадрат радиуса (для простоты радиусом я называю длину радиус-вектора, задающего положение случайной точки), то сначала выясним, как распределены радиусы. Так как круг заполнен равномерно, очевидно, что количество точек с радиусом r пропорционально длине окружности радиуса r. А длина окружности пропорциональна радиусу. Значит плотность распределения радиусов возрастает равномерно от центра окружности к её краям. А функция плотности имеет вид f(x) = 2x на интервале (0, 1). Коэффициент 2 для того, чтобы площадь фигуры под графиком равнялась единице. При возведении такой плотности в квадрат, она превращается в равномерную. Так как теоретически в данном случае для этого необходимо функцию плотности разделить на производную от функции преобразования (то есть от x²). А наглядно это происходит так:

Если аналогичное преобразование сделать для нормальной случайной величины, то функция плотности ее квадрата окажется похожей на гиперболу. А сложение двух квадратов нормальных случайных величин уже гораздо более сложный процесс, связанный с двойным интегрированием. И то, что в результате получится экспоненциальное распределение, лично мне тут остаётся проверить практическим методом или принять как аксиому. А кому интересно, предлагаю ознакомиться с темой поближе, почерпнув знаний из этих книжек:

• Вентцель Е.С. Теория вероятностей


• Кнут Д.Э. Искусство Программирования, том 2

В заключение приведу пример реализации генератора нормально распределенных случайных чисел на языке JavaScript:

function Gauss() {
        var ready = false;
        var second = 0.0;
        
        this.next = function(mean, dev) {
                mean = mean == undefined ? 0.0 : mean;
                dev = dev == undefined ? 1.0 : dev;
                
                if (this.ready) {
                        this.ready = false;
                        return this.second * dev + mean;
                }
                else {
                        var u, v, s;
                        do {
                                u = 2.0 * Math.random() - 1.0;
                                v = 2.0 * Math.random() - 1.0;
                                s = u * u + v * v;
                        } while (s > 1.0 || s == 0.0);
                        
                        var r = Math.sqrt(-2.0 * Math.log(s) / s);
                        this.second = r * u;
                        this.ready = true;
                        return r * v * dev + mean;
                }
        };
}

g = new Gauss(); // создаём объект
a = g.next(); // генерируем пару значений и получаем первое из них
b = g.next(); // получаем второе
c = g.next(); // снова генерируем пару значений и получаем первое из них

Параметры mean (математическое ожидание) и dev (среднеквадратическое отклонение) не обязательны. Обращаю ваше внимание на то, что логарифм натуральный и это мой первый пост.

This entry passed through the Full-Text RSS service — if this is your content and you're reading it on someone else's site, please read the FAQ at fivefilters.org/content-only/faq.php#publishers.