Как сэмплировать

Задача

Монте-Карло из предыдущих уроков предполагал, что сэмплировать из $p$ мы умеем. Обычно это не так: плотность известна с точностью до нормировки — например, $p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta)p(\theta)$, а знаменатель и есть тот интеграл, которого мы избегаем.

Значит, нужны методы, которым хватает возможности вычислять $p$ (или даже только $p$ с точностью до константы), а не сэмплировать из неё напрямую. Их три, и различаются они не столько механикой, сколько тем, что требуют и чем платят.

Виджет считает все три на одной двугорбой плотности на $[0,1]$. Двугорбой намеренно: на унимодальной цели сломанный сэмплер выглядит убедительно.

обратная CDFс отклонениемМетрополис

гистограмма и цель

сэмплов

2000

принято

2000 / 2000

доля принятых

100.0%

Обратная CDF

Если $U \sim \text{Uniform}(0,1)$, то $F^{-1}(U)$ распределено по $F$. Доказательство в одну строку: $P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x)$.

$$x = \htmlData{k=inv}{F^{-1}}\big(\htmlData{k=u}{u}\big), \qquad u \sim \text{Uniform}(0,1)$$

Метод точен, использует каждый сэмпл и не имеет настроек. Требование одно и оно жёсткое: нужна $F^{-1}$. Она есть у экспоненциального ($-\ln(1-u)/\lambda$), у Коши, у Лапласа — и её нет в замкнутой форме у гауссианы, и почти никогда нет в многомерном случае.

В виджете метод реализован табличным приближением: CDF считается на сетке один раз, дальше — двоичный поиск. Это тот же трюк, которым сэмплируют из категориального распределения по логитам, и в многомерном случае он не обобщается, потому что таблица растёт как $m^d$.

Отбор с отклонением

Возьмём огибающую $M q(x) \ge p(x)$. Сэмплируем $x \sim q$, $u \sim \text{Uniform}(0, Mq(x))$ и принимаем, если $u \le p(x)$. Геометрически: бросаем точки в прямоугольник и оставляем те, что попали под кривую.

Доля принятых равна ровно $1/M$, и это ахиллесова пята метода. Здесь $M = 2.9$ (пик нормированной плотности — $2.8039$, чуть меньше), и доля принятых выходит около $34\%$ — виджет показывает $33.1\%$ на двух тысячах черновых сэмплов.

Два способа испортить всё:

• огибающая не доминирует. Если $M q(x) < p(x)$ хоть где-то, горбы срезаются, и сэмплы приходят из неправильного распределения. Гистограмма при этом выглядит как распределение, и заметить подмену без независимой проверки нельзя;
• огибающая слишком высока. Метод остаётся корректным, но доля принятых падает, и вычисления уходят в мусор.

И основное ограничение: подобрать плотную огибающую в высокой размерности невозможно. Посчитать это можно точно. Огибаем $\mathcal{N}(0, I_d)$ распределением $\mathcal{N}(0, \sigma^2 I_d)$; отношение плотностей максимально в нуле и равно $M = \sigma^d$, значит доля принятых равна $\sigma^{-d}$:

$\sigma$	$d = 1$	$d = 10$	$d = 50$	$d = 100$
$1.1$	$0.91$	$0.39$	$8.5 \cdot 10^{-3}$	$7.3 \cdot 10^{-5}$
$1.5$	$0.67$	$1.7 \cdot 10^{-2}$	$1.6 \cdot 10^{-9}$	$2.5 \cdot 10^{-18}$
$2.0$	$0.50$	$9.8 \cdot 10^{-4}$	$8.9 \cdot 10^{-16}$	$7.9 \cdot 10^{-31}$

Обратите внимание, что даже огибающая, шире цели всего на $10\%$, к сотне измерений принимает один сэмпл из четырнадцати тысяч. А если промахнуться с масштабом вдвое — один из $10^{30}$. Метод не то чтобы медленный: он не работает.

Метрополис

Отказываемся и от CDF, и от огибающей. Идём случайным блужданием: из текущей точки предлагаем $x'$ симметрично, считаем

$\alpha = \min\!\left(1, \ \frac{p(x')}{p(x)}\right)$

и с вероятностью $\alpha$ переходим, иначе остаёмся на месте и записываем текущую точку снова.

Ключевое свойство — в дроби. Нормировка $p$ сокращается, поэтому Метрополис работает с ненормированной плотностью, и именно это делает его пригодным для байесовского вывода.

За это платят двумя вещами.

Сэмплы коррелированы. Отклонённый шаг дублирует предыдущую точку, а принятый уходит недалеко. Эффективное число независимых сэмплов много меньше их количества — та же величина $n_{\text{eff}}$, что была в уроке про importance sampling, только причина другая.

Есть настройка, и она нетривиальна. Подвигайте размер шага и смотрите не на долю принятых, а на гистограмму:

шаг	доля принятых	расстояние до цели ($n = 2000$)
$0.02$	$95\%$	$0.35$
$0.05$	$88\%$	$0.37$
$0.10$	$81\%$	$0.086$
$0.40$	$50\%$	$0.080$
$1.50$	$49\%$	$0.065$

Прочитайте первую строку внимательно: доля принятых $95\%$ — и худшая из всех гистограмма. Маленький шаг почти всегда принимается, потому что плотность в соседней точке почти такая же, но за две тысячи итераций цепь не успевает обойти носитель. Это и есть главное недоразумение вокруг MCMC: высокая доля принятых — симптом проблемы, а не признак здоровья.

Для сравнения, точный метод обратной CDF на тех же двух тысячах сэмплов даёт расстояние $0.053$ — это шум гистограммы, ниже которого не опустится никто. При шаге $0.4$ и двадцати тысячах сэмплов Метрополис доходит до $0.021$, то есть цепь действительно сходится — просто ей нужно больше итераций, чем независимому сэмплеру.

Что выбрать

	нужна CDF	нужна огибающая	работает с ненормированной $p$	размерность
обратная CDF	да	—	нет	1D практически
с отклонением	нет	да	да, если известна	низкая
Метрополис	нет	нет	да	высокая

Практика такая: обратной CDF сэмплируют из стандартных одномерных распределений и из категориального; отбор с отклонением используют внутри специализированных генераторов (в том числе гауссова — алгоритм зиккурата); всё остальное в байесовском выводе делает MCMC.

Современные варианты убирают именно случайность блуждания. HMC и NUTS ведут цепь по градиенту $\nabla \log p$, двигаясь вдоль плотности, а не наугад, и в высокой размерности выигрывают на порядки. Заметьте, что им нужен $\nabla \log p$ — та же score-функция, на которой стоит вся ветка про диффузионные модели. Ланжевеновская динамика, которую мы там встретим, — это буквально шаг градиентного подъёма по $\log p$ плюс гауссов шум, то есть Метрополис, которому подсказали направление.