Генератор нормальных значений: Бокс-Мюллер и полярный метод

Если выводишь распределение R² и понимаешь отбраковку точек квадрата, переходи к «Теперь сам».

Как думаешь, достаточно правильных среднего и дисперсии

Возьмём двенадцать независимых U(0,1), сложим и вычтем 6. Среднее равно нулю, потому что двенадцать половин сократились. Дисперсия равна 12·(1/12)=1. У стандартной нормали тоже среднее 0 и дисперсия 1. Значит ли это, что мы получили её точно?

Нет: сумма лежит между −6 и 6, а нормальный закон имеет ненулевые хвосты дальше любого конечного порога. Приближение может быть пригодно в одной задаче, но невозможность получить хвостовые значения важна в другой. Два момента не определяют распределение.

На что это похоже

Вместо выбора координат точки по отдельности выберем направление и расстояние от начала. Направление должно быть равномерным, но радиус нельзя взять равномерно от 0 до 1: дальние кольца имеют большую площадь, чем ближние.

У двумерной стандартной нормали плотность зависит только от расстояния. Круговая симметрия позволяет отделить угол от радиуса, после чего каждый выбирается простым одномерным правилом. Так двумерная задача становится удобнее двух отдельных обратных CDF, хотя изначально требовалась всего одна нормальная величина.

Механизм: обратная функция и полярные координаты

Для строго возрастающей непрерывной CDF F можно взять X=F⁻¹(U): событие X≤x равносильно U≤F(x), поэтому его вероятность F(x). Это общий способ генерации через обратную функцию распределения. Для нормальной Φ обратную функцию вычисляют численно; простого выражения через обычные элементарные функции нет.

У независимых стандартных X,Y совместная плотность равна exp(−(x²+y²)/2)/(2π). В круговом слое радиуса r площадь 2πr dr, поэтому плотность радиуса равна r exp(−r²/2) при r>0. Его хвост P(R>r)=exp(−r²/2). Угол независим от радиуса и равномерен на [0,2π): совместная плотность радиуса и угла факторизуется как r exp(−r²/2)·1/(2π).

Теперь зададим R=√(−2ln U₁). Для r≥0 событие R>r равносильно U₁<exp(−r²/2), поэтому вероятность имеет нужный хвост. Независимый U₂ задаёт θ=2πU₂. Возврат к x=Rcosθ,y=Rsinθ восстанавливает исходную совместную нормальную плотность, а значит две независимые координаты.

Правило

Преобразование Бокса-Мюллера использует два независимых равномерных источника внутри (0,1). Полярный метод берёт независимые a,b~U(−1,1), считает s=a²+b², принимает только 0<s<1 и возвращает a√(−2ln(s)/s), b√(−2ln(s)/s). После генерации стандартного Z нужная нормаль равна μ+σZ.

Две выдаваемые координаты используют общие R и θ, но всё же независимы: это следует из восстановленной совместной плотности. Общий источник промежуточных переменных сам по себе не доказывает зависимости итоговых величин.

Разбор: почему полярная отбраковка работает

Равномерная точка квадрата [−1,1]² после условия s<1 равномерна в единичном круге: плотность одинакова на принятой области, затем нормируется её площадью π. Вероятность принять попытку равна π/4. Среднее число попыток для пары равно 4/π, как у геометрического ожидания успеха.

Для точки круга P(s≤u)=πu/π=u при 0≤u≤1, поскольку круг радиуса √u имеет площадь πu. Значит s равномерно в [0,1]. Направление независимо от s. Нормированные a/√s,b/√s играют роли cosθ,sinθ, а √(−2ln s) даёт правильный нормальный радиус. Отсюда формула без вычисления тригонометрии.

import math
import random

rng = random.Random(95)


def box_muller():
    u = 1 - rng.random()
    angle = 2 * math.pi * rng.random()
    radius = math.sqrt(-2 * math.log(u))
    return radius * math.cos(angle), radius * math.sin(angle)


def polar():
    while True:
        a = 2 * rng.random() - 1
        b = 2 * rng.random() - 1
        s = a * a + b * b
        if 0 < s < 1:
            factor = math.sqrt(-2 * math.log(s) / s)
            return a * factor, b * factor


generators = [("Бокс-Мюллер", box_muller), ("полярный", polar)]
for label, generator in generators:
    samples = [x for _ in range(5000) for x in generator()]
    mean = sum(samples) / len(samples)
    variance = sum((x - mean) ** 2 for x in samples) / len(samples)
    print(label, "n", len(samples), "среднее", f"{mean:.3f}",
          "средний квадрат отклонения", f"{variance:.3f}")
print("12 равномерных минус 6: среднее 0, дисперсия 1, носитель [-6,6]")

Бокс-Мюллер n 10000 среднее 0.004 средний квадрат отклонения 1.038
полярный n 10000 среднее 0.023 средний квадрат отклонения 0.989
12 равномерных минус 6: среднее 0, дисперсия 1, носитель [-6,6]

Строки показывают конечные выборки по 10 000 значений. Их средние и средние квадраты отклонений близки к 0 и 1, но не обязаны совпасть. Это проверка реализации на выборке, а не доказательство распределения; доказательство дано через совместную плотность. Random(95) фиксирует последовательность для повторяемого запуска.

Разбор: границы машинного источника

В коде rng.random() выдаёт число от 0 включительно до 1 не включительно. В Боксе-Мюллере берём 1−rng.random(), поэтому аргумент логарифма положителен. В идеальной непрерывной модели включение одного конца не меняет закон. На машине множество значений конечно, и оно даёт приближение, а не буквальную бесконечную точность нормальной случайной величины.

В полярном методе условие 0<s<1 исключает центр, где появились бы ln(0) и деление на ноль, и точки вне круга. Отбраковка нужна не только ради защиты от ошибки функции: без неё s имело бы другой закон, и распределение результата изменилось бы. Заменить неудачный s на 1 также нельзя: так появились бы лишние нулевые пары вместо независимых нормальных значений.

Если нужны N(μ,σ²), умножают результат на σ, а не на σ². Для многомерной модели независимые стандартные значения сначала объединяют в Z, затем вычисляют AZ+μ. Случайно повторить один и тот же Z во всех координатах означало бы создать совсем другую, вырожденную ковариацию.

Что ломается без этого

Путают равномерность радиуса и равномерность площади. Равномерный радиус концентрирует слишком большую долю точек около центра: площадь внутреннего круга растёт как квадрат радиуса. Поэтому в отбраковке равномерна величина s=r², а не r при равномерной точке круга.

Также нельзя подтверждать нормальность только средним и дисперсией. Пример суммы двенадцати равномерных проходит такую проверку, но не способен выдавать значения за пределами [−6,6]. Нужны проверка механизма и анализ области возможных результатов, особенно если задача зависит от редких больших отклонений.

Зачем это в работе

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

Теперь сам

Какова вероятность принять одну попытку полярного метода? Как получить N(5,9) из стандартного Z? Верна ли для суммы двенадцати равномерных вероятность получить число 7?