Метод отклонения (rejection sampling)

Метод обратной функции элегантен, но капризен: он требует, чтобы функцию распределения можно было «развернуть» формулой. А что делать, когда CDF не обращается — или когда вы знаете лишь форму плотности, но не её интеграл? На помощь приходит грубоватый, но удивительно универсальный приём — метод отклонения. Его идея настолько проста, что её можно объяснить, бросая дротики в мишень.

Метод отклонения (rejection sampling) — способ генерации выборки из распределения с плотностью f(x): случайные точки разбрасывают по прямоугольнику, накрывающему график плотности, и оставляют лишь те, что попали под кривую; их x-координаты и образуют нужную выборку.

Метафора «дротиков» здесь буквальна. Представьте, что график плотности — это холм, нарисованный внутри прямоугольной рамки. Вы вслепую кидаете дротики, равномерно покрывая всю рамку. Дротики, воткнувшиеся под линией холма, мы засчитываем; угодившие над ней — выбрасываем. Там, где холм высокий, под кривой много места — туда попадает больше засчитанных дротиков; где холм низкий — мало. В итоге горизонтальные координаты принятых дротиков воспроизводят форму плотности, какой бы причудливой она ни была.

Когда обратная функция бессильна

Главное достоинство метода отклонения — ему не нужна ни CDF, ни её обратная функция. Достаточно уметь вычислять саму плотность f(x) в любой точке. Это огромное преимущество: многие плотности (особенно в байесовской статистике и физике) заданы сложными формулами, интеграл от которых в элементарных функциях не берётся, а значит, метод обратной функции к ним неприменим. Метод отклонения же работает «в лоб», требуя лишь умения сказать «вот здесь кривая на такой-то высоте».

Разбор на плотности f(x) = 2x

Возьмём конкретную целевую плотность: f(x) = 2x на отрезке [0, 1]. Это корректная плотность — площадь под ней равна 1 (треугольник с катетами 1 и 2 даёт площадь ½·1·2 = 1). Распределение «перекошено» вправо: большие x вероятнее малых, ведь высота кривой растёт линейно.

Построим обрамляющий прямоугольник. По горизонтали берём весь носитель [0, 1]. По вертикали нужно накрыть максимум плотности: f(x) = 2x на [0, 1] достигает наибольшего значения 2 в точке x = 1. Значит, прямоугольник — это [0, 1] × [0, 2]. Алгоритм:

1. Сгенерировать x равномерно на [0, 1] и u равномерно на [0, 2] — это случайная точка в прямоугольнике.
2. Если u ≤ f(x), то есть u ≤ 2x — точка под кривой, принимаем x.
3. Иначе отбрасываем точку и пробуем снова.

import random, statistics
random.seed(4)
samples = []
tries = 0
while len(samples) < 10000:
    tries += 1
    x = random.random()
    u = random.random() * 2
    if u <= 2 * x:
        samples.append(x)
print(f"Принято {len(samples)} из {tries}")
print(f"Доля принятых: {len(samples) / tries:.2%}")
print(f"Среднее (теория 2/3 = 0.6667): {statistics.mean(samples):.4f}")

Вывод:

Принято 10000 из 19996
Доля принятых: 50.01%
Среднее (теория 2/3 = 0.6667): 0.6678

Разберём результаты. Чтобы набрать 10000 принятых выборок, понадобилось 19996 попыток — почти вдвое больше. Доля принятых составила 50.01%, то есть ровно половину дротиков мы выбросили. И главное: среднее принятых x получилось 0.6678 — практически совпало с теоретическим средним этого распределения 2/3 ≈ 0.6667. Метод честно воспроизвёл перекос вправо: раз большие x под кривой имеют больше «места», среднее сместилось от центра 0.5 к 0.667.

Как работает под капотом

Почему x-координаты принятых точек распределены именно по f(x)? Ключ — в том, что мы разбрасываем точки равномерно по площади прямоугольника. Рассмотрим узкую вертикальную полоску шириной dx около какого-то x. Высота «принимаемой» зоны в этой полоске равна f(x) (от 0 до кривой), значит, доля принятых точек, попавших в полоску, пропорциональна f(x). А «пропорционально плотности» — это ровно и есть определение выборки из распределения f. Никакой магии: равномерность по площади плюс отсев под кривой автоматически дают нужную форму.

Теперь об эффективности — почему приняли именно около 50%. Доля принятых точек равна отношению площадей: площадь под кривой к площади прямоугольника. Площадь под нашей плотностью равна 1 (это свойство любой плотности). Площадь прямоугольника [0, 1] × [0, 2] равна 1 * 2 = 2. Отношение 1/2 = 50% — ровно то, что показал эксперимент. Отсюда практический вывод: чем плотнее прямоугольник облегает кривую, тем меньше точек пропадает зря. Если бы максимум плотности был не 2, а 200, прямоугольник раздулся бы, и эффективность рухнула бы до 0.5% — пришлось бы делать в сотни раз больше попыток ради той же выборки. Это и есть главная слабость метода: в «худо подобранной» рамке он расточителен.

Частые ошибки

• Сделать прямоугольник ниже максимума плотности. Если верхняя граница по u меньше пика f(x), часть кривой окажется выше рамки и «срежется» — выборка получится искажённой. Высоту прямоугольника берут не меньше максимума плотности.
• Раздувать прямоугольник без нужды. Слишком высокая или широкая рамка ничего не ломает, но резко роняет эффективность: большинство точек отбрасывается, и метод работает мучительно долго. Рамку стараются подогнать вплотную к кривой.
• Принимать точку по неверному сравнению. Условие принятия — u ≤ f(x). Перепутать знак (u ≥ f(x)) или сравнить не с той величиной — значит получить «зеркальное» или вовсе бессмысленное распределение.
• Забыть, что число попыток заранее неизвестно. Сколько итераций уйдёт на нужное число принятых выборок — случайная величина. Цикл должен крутиться до набора выборки (while len(samples) < N), а не фиксированное число раз, иначе наберётся меньше, чем нужно.
• Считать отброшенные точки «потерей качества». Отбрасывание — не дефект, а сам механизм метода. Принятые точки абсолютно корректны; выброшенные просто не нужны. Расточительность бьёт по скорости, но не по правильности.

Итоги

• Метод отклонения генерирует выборку из плотности f(x), не требуя ни CDF, ни обратной функции — достаточно уметь вычислять саму плотность.
• Точки бросают равномерно в прямоугольник, накрывающий график плотности, и принимают те, для которых u ≤ f(x).
• Для f(x) = 2x на [0, 1] рамка — [0, 1] × [0, 2], среднее принятых сходится к 2/3.
• Эффективность = (площадь под кривой) / (площадь прямоугольника); для нашего примера это 1/2, поэтому принимается около 50% точек.
• Чем плотнее рамка облегает кривую, тем меньше точек пропадает; плохо подобранная рамка делает метод расточительным, но не неправильным.