Методы Монте-Карло

Закон больших чисел превратил любое математическое ожидание в предел средних: берём сэмплы, усредняем, и ошибка убывает как $1/\sqrt n$. Чтобы этот рецепт стал методом, нужно ответить ещё на два вопроса. Как вообще получить сэмплы из заданного распределения, если компьютер выдаёт только равномерно распределённые случайные числа? И что делать, когда $1/\sqrt n$ — это слишком медленно, как в случае редких событий? Эта заметка отвечает на оба вопроса: метод обратной функции, преобразования и выборка с отклонением дают сэмплы, а выборка по значимости и снижение дисперсии делают каждый сэмпл ценнее.

Обозначения. Цель — величина $\mu = E_p[g(X)]$ для плотности $p$; $\hat\mu_n = \frac1n \sum g(X_i)$ при $X_i \sim p$. Предлагающее распределение (proposal) $q$ — плотность, из которой легко сэмплировать.

Простой метод Монте-Карло и его ошибка#

Если $X_i$ независимы и имеют плотность $p$, а $E_p\lvert g(X) \rvert < \infty$, то оценка $\hat\mu_n$ несмещённая и сходится к $\mu$. Если к тому же $\sigma^2 = \operatorname{Var}_p(g(X))$ конечна, то по центральной предельной теореме она приближённо распределена как $N(\mu, \sigma^2/n)$. Оценку методом Монте-Карло всегда нужно приводить вместе с интервалом погрешности $\hat\mu_n \pm 1.96\, \hat\sigma/\sqrt n$, вычисленным по тем же сэмплам. Размерность $X$ не входит в ошибку напрямую: ошибка определяется только $\sigma$ и $n$, поэтому для интегралов высокой размерности и выбирают метод Монте-Карло. Без конечной дисперсии интервал погрешности бессмыслен, даже если оценка сходится: для $g(X)$ с тяжёлыми хвостами и для весов значимости с бесконечной дисперсией ниже нужны другие диагностики или устойчивые оценки, например медиана средних (median of means). Всё, что следует ниже, либо порождает сэмплы, либо уменьшает $\sigma$.

Как получить сэмплы#

Генератор случайных чисел выдаёт независимые числа с распределением $U(0, 1)$. Три приёма превращают их в сэмплы из других распределений.

• Метод обратной функции (inverse transform sampling). $F^{-1}(U)$ имеет функцию распределения $F$, как показано в заметке о функциях от случайных величин: $-\ln(1 - U)$ имеет показательное (экспоненциальное) распределение, а поиск по накопленным суммам даёт сэмплы из любого дискретного распределения.
• Преобразования. Когда у $F^{-1}$ нет формулы, она может найтись у другого отображения. Преобразование Бокса–Мюллера (Box–Muller transform) переводит две независимые равномерные величины в стандартную нормальную: $Z = \sqrt{-2\ln(1 - U_1)} \cos(2\pi U_2)$. Оно работает потому, что у пары $(Z_1, Z_2)$ независимых нормальных величин плотность инвариантна относительно поворотов, поэтому угол этой точки распределён равномерно и не зависит от радиуса, а квадрат радиуса имеет показательное распределение со средним 2.
• Выборка с отклонением (rejection sampling). Если плотность известна только с точностью до константы или для неё нет никакой удобной формулы сэмплирования, сэмплы предлагают, а затем принимают или отклоняют.

Выборка с отклонением#

Выберем предлагающее распределение $q$, из которого легко сэмплировать, и константу $M$, такую что $p(x) \le M q(x)$ всюду. Сэмплируем $X \sim q$ и $U \sim U(0, 1)$ и принимаем $X$, если $U M q(X) \le p(X)$; иначе отбрасываем его и пробуем снова. Точка $(X, U M q(X))$ распределена равномерно под огибающей (envelope) $M q$, поэтому принятые точки распределены равномерно под графиком $p$, а их $x$-координаты имеют плотность $p$:

Рисунок

Выборка с отклонением

Методу нужна $p$ только с точностью до постоянного множителя, поскольку константу можно включить в $M$. Его цену задаёт доля принятых (acceptance rate), для нормированной плотности равная $1/M$: в среднем на один принятый сэмпл приходится $M$ предложений. В примере это $2.46$, но с ростом размерности это число, как правило, растёт экспоненциально. При сэмплировании единичного шара с отклонением из описанного вокруг него куба принимается 79% предложений в двумерном случае, 0.25% в десятимерном и $2.5 \cdot 10^{-8}$ в двадцатимерном: почти весь объём многомерного куба сосредоточен в его углах — то же явление заметка о многомерном нормальном распределении обнаружила для оболочки гауссианы. Выборка с отклонением — инструмент для малых размерностей и для внутренних циклов других сэмплеров.

Выборка по значимости#

Простой метод Монте-Карло тратит сэмплы там, где велика $p$, а это может быть далеко от области, где важна $g$. Выборка по значимости (importance sampling) вместо этого сэмплирует из предлагающего распределения $q$ и исправляет результат весами значимости (importance weights):

$$E_p[g(X)] = \int g(x) \frac{p(x)}{q(x)}\, q(x)\, dx = E_q\big[g(X)\, w(X)\big], \qquad w(x) = \frac{p(x)}{q(x)} ,$$

что верно, если $q > 0$ всюду, где $g p \ne 0$. Оценка $\frac1n \sum g(X_i) w(X_i)$ при $X_i \sim q$ несмещённая, а её дисперсия зависит от того, насколько хорошо $q$ соответствует $\lvert g \rvert p$; оптимальное предлагающее распределение пропорционально этому произведению.

Сильнее всего выигрыш для редких событий. Чтобы оценить $P(Z > 4) = 3.17 \cdot 10^{-5}$ для стандартной нормальной $Z$, миллион простых сэмплов дают всего 28 событий при ожидаемых $31.7$: относительная стандартная ошибка около 18%, а фактическая ошибка в этом запуске — 12%. Сэмплирование из $q = N(4, 1)$ помещает в событие половину сэмплов, и каждый из них получает вес $p/q = e^{8 - 4x}$:

Рисунок

Выборка по значимости для редкого события

Десять тысяч сэмплов из $q$ дают оценку $3.27 \cdot 10^{-5}$ с относительной стандартной ошибкой 2% — в сто раз меньше сэмплов при ошибке в восемь раз меньше; фактическая ошибка этой оценки, известная здесь потому, что известно точное значение, равна 3.3%. Для других вопросов те же веса бесполезны: их эффективный размер выборки (effective sample size, ESS) $(\sum w_i)^2 / \sum w_i^2$ равен $5.1$ из 10 000, так что оценивать с ними, скажем, среднее $Z$ было бы безнадёжно. Предлагающее распределение хорошо лишь для конкретной подынтегральной функции, а предлагающее распределение с более лёгкими хвостами, чем у $p$, может давать бесконечную дисперсию.

Когда $p$ известна только с точностью до константы, как апостериорное распределение, веса нормируют на их сумму (самонормированная выборка по значимости, self-normalized importance sampling); тогда оценка немного смещена, но состоятельна. Веса значимости встречаются во всём машинном обучении: оценка политики по данным другой политики (off-policy evaluation) в обучении с подкреплением перевзвешивает записанные действия отношением новой политики к старой, а нижняя оценка с весами значимости (importance-weighted bound) для вариационных автоэнкодеров усредняет веса значимости нескольких сэмплов.

Снижение дисперсии#

Даже при правильном распределении дисперсию оценки часто можно уменьшить, используя то, что уже известно.

• Контрольные переменные (control variates). Если величина $h(X)$ с известным средним коррелирует с $g(X)$, будем оценивать вместо исходного математического ожидания $E[g(X) + c\,(h(X) - E h(X))]$ с $c = -\operatorname{Cov}(g, h)/\operatorname{Var}(h)$. Чтобы вычислить $E[e^U] = e - 1$ для $U \sim U(0, 1)$, контрольная переменная $h(U) = U$ со средним $1/2$ снижает дисперсию на один сэмпл с $0.2446$ до $0.0040$, в 60 раз, без дополнительных затрат. В обучении с подкреплением методами градиента политики (policy gradient) вычитание базовой линии (baseline) из вознаграждения — это контрольная переменная для оценки градиента, и это одна из главных причин, по которым такие методы работают на практике.
• Антитетические пары (antithetic). Используем $U$ и $1 - U$ вместе: для монотонной $g$ два значения отрицательно коррелированы, и их среднее имеет меньшую дисперсию, чем среднее двух независимых сэмплов.
• Стратификация и квазимонте-карло (stratification, quasi-Monte Carlo). Сэмплы распределяют по пространству равномерно, а не случайно; в малых размерностях последовательности с малой дискрепансией (low-discrepancy sequences) дают ошибку, близкую к $1/n$, а не к $1/\sqrt n$.

Оценка градиентов в машинном обучении — та же задача. Градиент математического ожидания $E_{q_\theta}[f(Z)]$ можно оценить через градиент логарифма плотности — это оценка REINFORCE (score-function estimator), несмещённая, но с большой дисперсией, — или с помощью трюка репараметризации (reparameterization trick) из вариационного автоэнкодера, который записывает $Z = \mu_\theta + \sigma_\theta \varepsilon$ и обычно даёт гораздо меньшую дисперсию.

Что изменилось#

	Простой метод Монте-Карло	Методы Монте-Карло
Сэмплы из $p$	предполагаются доступными	метод обратной функции, преобразования, выборка с отклонением
Ошибка	$\sigma/\sqrt n$, приводится вместе с оценкой	та же скорость, меньшая $\sigma$
Редкие события	большинство сэмплов тратится впустую	выборка по значимости из предлагающего распределения около события
Известная структура	не используется	контрольные переменные, антитетические пары, стратификация
Высокая размерность	скорость не зависит от $d$	но доля принятых и веса значимости ухудшаются с ростом $d$

Ограничения#

• В высокой размерности трудно найти хорошее предлагающее распределение. И выборке с отклонением, и выборке по значимости нужно $q$, близкое к $p$; в многомерном случае даже небольшое несоответствие по каждой координате накапливается, и доля принятых и эффективный размер выборки падают экспоненциально.
• Выборка по значимости может отказать незаметно. Предлагающее распределение с более лёгкими хвостами, чем у целевого, может давать оценки с бесконечной дисперсией: они долго выглядят устойчивыми, а потом скачут.
• Точные сэмплы из сложных распределений доступны редко. Апостериорным распределениям реальных моделей, энергетическим моделям (energy-based models) и распределениям на структурах нужен другой подход — ему посвящена следующая заметка.

Что нужно запомнить#

• Оценку методом Монте-Карло приводят вместе с интервалом погрешности $\pm 1.96\,\hat\sigma/\sqrt n$, верным для независимых сэмплов, когда дисперсия $g(X)$ конечна; скорость не зависит от размерности.
• Сэмплы получают методом обратной функции ($F^{-1}(U)$), преобразованиями (Бокса–Мюллера) или выборкой с отклонением (принять $X \sim q$ с вероятностью $p(X)/(Mq(X))$, доля принятых $1/M$).
• Выборке с отклонением нужна $p$ лишь с точностью до константы, но в высокой размерности доля принятых обрушивается.
• Выборка по значимости: $E_p[g] = E_q[g\, p/q]$; она может сделать редкие события дешёвыми; нужно проверять эффективный размер выборки $(\sum w)^2/\sum w^2$.
• Контрольные переменные вычитают коррелированную величину с известным средним; базовые линии в методах градиента политики — это контрольные переменные.

Код#

Выборка с отклонением, преобразование Бокса–Мюллера, редкое событие через выборку по значимости и контрольная переменная:

import numpy as np, math
from scipy import stats

rng = np.random.default_rng(0)

# rejection sampling from Beta(2, 5) with a uniform proposal and M = max p = p(0.2)
p = stats.beta(2, 5).pdf
M = p(0.2)
x, u = rng.random(100_000), rng.random(100_000)
sample = x[u * M <= p(x)]
print(round(float(M), 4), round(len(sample) / 100_000, 3), round(float(sample.mean()), 3), round(2 / 7, 3))   # 2.4576 0.408 0.286 0.286

# the share of a cube taken by the inscribed ball: the acceptance rate of rejection sampling in d dimensions
print([f'{math.pi ** (d / 2) / math.gamma(d / 2 + 1) / 2 ** d:.2g}' for d in (2, 5, 10, 20)])   # ['0.79', '0.16', '0.0025', '2.5e-08']

# Box-Muller: two uniforms give a standard normal
u1, u2 = rng.random((2, 500_000))
z = np.sqrt(-2 * np.log(1 - u1)) * np.cos(2 * np.pi * u2)
print(round(float(z.mean()), 3), round(float(z.std()), 3))           # 0.0 1.003

# a rare event, P(Z > 4): a million plain samples against 10,000 importance samples from N(4, 1)
true = stats.norm.sf(4)
plain = rng.normal(size=1_000_000) > 4
y = rng.normal(4, 1, 10_000)
w = np.exp(8 - 4 * y)                                                # p(y) / q(y)
terms = (y > 4) * w
print(f'{true:.3g}', int(plain.sum()), f'{terms.mean():.3g}', round(float(terms.std() / math.sqrt(1e4) / true), 3))   # 3.17e-05 28 3.27e-05 0.022
print(round(float(w.sum() ** 2 / (w ** 2).sum()), 1))               # 5.1: the effective sample size of these weights

# a control variate for E[e^U] = e - 1: the variance per sample without and with it
U = rng.random(10_000)
Y = np.exp(U)
c = -np.cov(Y, U)[0, 1] / U.var(ddof=1)
print(round(float(Y.mean()), 4), round(float(Y.var()), 4), round(float((Y + c * (U - 0.5)).var()), 4), round(math.e - 1, 4))   # 1.7181 0.2446 0.004 1.7183

Дальше#

Выборке с отклонением и выборке по значимости нужно предлагающее распределение, близкое к целевому, а в высокой размерности такого обычно нет. Последняя заметка курса вместо этого строит цепь Маркова, стационарное распределение которой и есть целевое, используя только ненормированную плотность: алгоритм Метрополиса–Гастингса, сэмплирование по Гиббсу и динамику Ланжевена — мост к диффузионным моделям.

Литература#

• C. P. Robert, G. Casella. Monte Carlo Statistical Methods, 2nd ed. Springer, 2004. Глава 2 (генерация случайных величин, метод отклонения), глава 3 (интегрирование методом Монте-Карло, выборка по значимости), глава 4 (снижение дисперсии).
• A. B. Owen. Monte Carlo Theory, Methods and Examples. 2013. Главы 4, 8 и 9 (неравномерно распределённые случайные величины, снижение дисперсии, выборка по значимости).
• G. E. P. Box, M. E. Muller. A Note on the Generation of Random Normal Deviates. The Annals of Mathematical Statistics, 29(2), 1958.
• J. von Neumann. Various Techniques Used in Connection with Random Digits. National Bureau of Standards Applied Mathematics Series, 12, 1951. Метод отклонения.
• R. J. Williams. Simple Statistical Gradient-Following Algorithms for Connectionist Reinforcement Learning. Machine Learning, 8, 1992. REINFORCE и базовые линии.