Урок 67. MCMC и расшифровка текста

MCMC не умеет напрямую вытянуть число из сложного распределения. Вместо этого метод строит зависимую прогулку, которая достаточно долго живёт в нужных областях. Из локального правила перехода возникает глобальная статистика.

Когда плотность известна только с точностью до множителя

В байесовской задаче апостериорная плотность параметра $\theta$ имеет вид

p(\theta\mid D)= \frac{p(D\mid\theta)p(\theta)} {\int p(D\mid u)p(u)\,du}.

Числитель вычислить часто можно, а многомерный интеграл в знаменателе — нет. Обозначим ненормированную плотность

\widetilde\pi(\theta)=p(D\mid\theta)p(\theta), \qquad \pi(\theta)=\frac{\widetilde\pi(\theta)}{Z}.

Для ожидания $I=\mathbb E_\pi f(\theta)$ нужны выборки из $\pi$ , но стандартный генератор требует нормированного распределения. MCMC обходит круг: строит цепь Маркова с заданным стационарным распределением $\pi$ и оценивает

\widehat I_N=\frac1N\sum_{t=1}^N f(\theta_t).

Выборки $\theta_t$ зависимы. Цена отказа от неизвестной константы $Z$ — необходимость анализировать перемешивание цепи. Базовые свойства стационарности и переходов были введены в цепях Маркова, а теперь они становятся вычислительным инструментом.

Двухвершинная плотность, путь цепи Метрополиса и гистограмма её состояний — Рис. 67.1. Целевая плотность и зависимая траектория
Вверху показана ненормированная двухмодальная плотность $\widetilde\pi(x)$ . В середине — первые 400 состояний цепи с редкими переходами между модами. Внизу гистограмма после прогрева сопоставлена с нормированной $\pi$ ; совпадение маргинальных частот ещё не гарантирует хорошего перемешивания.

Баланс потоков

Пусть $K(x,dy)$ — переходное правило. Достаточное условие стационарности называется детальным балансом:

\pi(x)K(x,dy)=\pi(y)K(y,dx).

Оно говорит, что в равновесии поток вероятности из $x$ в $y$ компенсируется обратным потоком. Если просуммировать по всем $x$ , входящий поток в $y$ станет равен $\pi(y)$ . Детальный баланс сильнее необходимого: существуют несимметричные цепи со стационарным $\pi$ , но обратимость удобно доказывать.

В дискретном случае это равенство потоков по каждой паре состояний. В непрерывном — равенство мер; плотности можно читать аналогично, помня про элементы объёма.

Алгоритм Метрополиса–Хастингса

Находясь в $x$ , предложим новую точку $y\sim q(y\mid x)$ . Примем её с вероятностью

a(x,y)=\min\left( 1,\; \frac{\widetilde\pi(y)q(x\mid y)} {\widetilde\pi(x)q(y\mid x)} \right).

Если предложение отклонено, следующее состояние снова равно $x$ . Повторы обязательны: без них стационарное распределение изменится. Неизвестная константа $Z$ сокращается в отношении.

Для симметричного случайного предложения $q(y\mid x)=q(x\mid y)$ остаётся

a(x,y)=\min\left(1,\frac{\widetilde\pi(y)}{\widetilde\pi(x)}\right).

Движение в область большей плотности принимается всегда, а вниз по плотности — иногда. Эти редкие спуски позволяют выйти из локального пика. Формула похожа на критерий в имитации отжига, однако цель MCMC — сохранить распределение, а не найти единственный максимум.

Масштаб предложения и цена корреляции

Рассмотрим random-walk Metropolis:

y=x+\sigma\varepsilon,\qquad \varepsilon\sim\mathcal N(0,I).

При слишком малом $\sigma$ почти всё принимается, но цепь ползёт микроскопическими шагами. При слишком большом $\sigma$ предложения падают в области ничтожной плотности и отклоняются; цепь снова почти стоит. Между крайностями есть рабочий масштаб, зависящий от размерности и геометрии $\pi$ .

Автокорреляция функции $f(\theta_t)$ на лаге $k$ :

\rho_k=\operatorname{Corr}\bigl(f(\theta_t),f(\theta_{t+k})\bigr).

Интегрированное время автокорреляции

\tau_{\mathrm{int}}=1+2\sum_{k\ge1}\rho_k

уменьшает эффективный размер выборки:

N_{\mathrm{eff}}\approx\frac{N}{\tau_{\mathrm{int}}}.

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

Трассы, доли принятия и автокорреляции цепей с малым, подходящим и большим шагом — Рис. 67.2. Три масштаба предложения
Столбцы соответствуют $\sigma=0{,}05$ , $0{,}8$ и $6$ . В каждой колонке показаны трасса, доля принятых предложений и $\rho_k$ . Малый шаг даёт высокую acceptance rate, большой — длинные повторы; средний быстрее разрушает корреляцию.

Лаборатория цепи

Выберите двухвершинную целевую плотность и запустите несколько цепей из левой моды, правой моды и промежутка. Меняйте ширину предложения. Записывайте долю принятия, число переходов между модами, среднее по каждой цепи и эффективный размер выборки.

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

Прогрев, несколько цепей и диагностика

Начальный участок, пока цепь забывает старт, называют warm-up или burn-in. Просто удалить первые 10 процентов недостаточно: нужная длина зависит от переходов. Лучше запускать несколько цепей из разнесённых точек и сравнивать их.

Диагностика $\widehat R$ сопоставляет разброс между цепями и внутри них. Значение около единицы — необходимый, но не достаточный признак. Цепи могут вместе застрять в одной моде. Нужны трассы, эффективный размер выборки, автокорреляции и проверки на искусственных данных, где ответ известен.

Thinning, то есть сохранение каждого $k$ -го состояния, обычно не создаёт информацию. Оно уменьшает файл, но выбрасывает вычисленные точки. Для оценки среднего выгоднее учитывать все состояния и корректно оценивать автокорреляцию.

Энергия, температура и многомодальность

Плотность часто записывают как

\pi(x)\propto e^{-E(x)/T}.

Энергия $E$ мала в вероятных состояниях, температура $T$ сглаживает рельеф. При низкой температуре пики узки, переходы между ними редки. Tempering-методы запускают цепи при разных температурах и обменивают состояния: горячая цепь пересекает барьеры, холодная сохраняет нужную цель.

Связь с физикой не декоративна. Формула Больцмана дала язык и для статистической механики, и для вероятностных алгоритмов. Позже похожая экспоненциальная нормировка появится в attention как softmax, хотя смысл переменных будет другим.

Реальный пример: байесовская логистическая модель

Возьмём открытый датасет UCI Heart Disease. Цель $y_i\in\{0,1\}$ обозначает наличие диагноза, признаки $x_i$ включают возраст, давление и лабораторные показатели. После стандартизации зададим

\Pr(y_i=1\mid x_i,\beta)=\sigma(x_i^\top\beta), \qquad \beta_j\sim\mathcal N(0,s^2).

Логарифм ненормированного posterior:

\log\widetilde\pi(\beta) =\sum_i\left[ y_i\log\sigma(x_i^\top\beta) +(1-y_i)\log(1-\sigma(x_i^\top\beta)) \right] -\frac{\|\beta\|^2}{2s^2}.

Random-walk Metropolis можно реализовать в нескольких строках, но коррелированные признаки создают вытянутую геометрию posterior. Одинаковый шаг по всем координатам работает плохо. Масштабирование признаков и предложение с ковариацией, приближённой к posterior, резко улучшают перемешивание.

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

Контуры апостериорной плотности двух коэффициентов и траектории изотропного и адаптированного MCMC — Рис. 67.3. Коррелированный posterior логистической модели
Контуры показывают вытянутую совместную плотность двух коэффициентов. Красная цепь с изотропным предложением движется поперёк узкой долины и часто отклоняется; синяя использует ориентированную ковариацию. Рядом указаны $N_{\mathrm{eff}}$ для обоих коэффициентов при одинаковом числе вычислений.

Мини-исследование: проверить MCMC на задаче с известным ответом

Перед сложным posterior полезно устроить калибровочный стенд. Возьмите двумерную normal target с известными средним $(0,0)$ и covariance

\Sigma= \begin{pmatrix} 1&0{,}95\\ 0{,}95&1 \end{pmatrix}.

Запустите random-walk Metropolis с изотропным предложением и с предложением, covariance которого пропорциональна $\Sigma$ . Для каждого варианта используйте четыре разнесённых старта и одинаковое число вычислений $\widetilde\pi$ . Сопоставьте acceptance rate, ошибку среднего, covariance и $N_{\mathrm{eff}}$ по направлениям $(1,1)$ и $(1,-1)$ .

Узкое направление $(1,-1)$ и длинное $(1,1)$ имеют разные масштабы. Изотропный шаг, подходящий первому, медленно движется вдоль второго; шаг для длинного направления часто вылетает поперёк. Ориентированное предложение учитывает геометрию. Если диагностика не обнаруживает эту разницу на известной цели, ей рано доверять в реальной модели.

Повторите опыт после линейного whitening-преобразования $u=\Sigma^{-1/2}x$ . В $u$ -координатах target сферична, и простой шаг должен работать лучше. Это вычислительная версия стандартизации признаков из линейных моделей: смена координат не меняет вероятностный вопрос, но меняет трудность прогулки.

Для каждой оценки среднего постройте Monte Carlo standard error

\operatorname{MCSE}(\bar f) \approx\sqrt{\widehat{\operatorname{Var}}(f)/N_{\mathrm{eff}}}.

Проверьте coverage: в каком проценте независимых запусков интервал $\bar f\pm1{,}96\,\mathrm{MCSE}$ содержит известное истинное значение. Диагностика становится проверяемой процедурой, а не набором красивых trace plots.

Три слоя успешной цепи

Успешная MCMC-работа состоит из трёх слоёв. Математика гарантирует нужное стационарное распределение; вычислительная диагностика проверяет, исследовала ли цепь его существенные области; предметная проверка спрашивает, осмысленна ли сама вероятностная модель. Доля принятия относится только ко второму слою и не может заменить остальные.

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

Задачи

3баллов

На состояниях $(A,B,C)$ ненормированные веса цели равны $(1,2,7)$ . Из текущего состояния предложение выбирает каждое из двух других с вероятностью $1/2$ , а принятие задаёт $a(x,y)=\min\{1,\widetilde\pi(y)/\widetilde\pi(x)\}$ . Выпишите $a(x,y)$ для всех шести направлений и строко-стохастическую матрицу перехода, включая вероятности остаться. Начальное состояние и seed не нужны; зачёт требует точные дроби и проверку суммы каждой строки.

4баллов

Пусть $q(y\mid x)=q(x\mid y)$ , $\pi(x)=\widetilde\pi(x)/Z$ и $a(x,y)=\min\{1,\widetilde\pi(y)/\widetilde\pi(x)\}$ . Докажите равенство $\pi(x)q(y\mid x)a(x,y)=\pi(y)q(x\mid y)a(y,x)$ отдельно для случаев $\widetilde\pi(y)\ge\widetilde\pi(x)$ и $\widetilde\pi(y)<\widetilde\pi(x)$ . Затем просуммируйте по $x$ и получите стационарность $\pi$ . Случайности, начального состояния и единиц измерения нет; в полном решении должно быть явно показано сокращение неизвестного $Z$ .

5баллов

Для цели $\pi(x)=\mathcal N(0,1)$ запустите четыре random-walk Metropolis цепи из общего старта $x_0=8$ с предложением $y=x+\sigma\varepsilon$ , $\varepsilon\sim\mathcal N(0,1)$ и $\sigma\in\{0{,}05;0{,}5;2;10\}$ . В каждом запуске с PCG64 seed 6703+j, где $j=0,1,2,3$ в порядке значений $\sigma$ , сделайте 60 000 итераций, отбросьте первые 10 000 как warm-up и анализируйте оставшиеся $N=50\,000$ . Оцените $\tau_{\mathrm{int}}=1+2\sum_{k=1}^{K}\widehat\rho_k$ методом начальной положительной последовательности Гейера: суммируйте пары $\widehat\rho_{2m-1}+\widehat\rho_{2m}$ до первой неположительной пары; $N_{\mathrm{eff}}=N/\tau_{\mathrm{int}}$ . Зачёт требует acceptance rate, среднее, несмещённую дисперсию, $\widehat\rho_1$ , $\tau_{\mathrm{int}}$ , $N_{\mathrm{eff}}$ и воспроизводимые трассы.

6баллов

Цель равна $0{,}5\mathcal N(-4,1)+0{,}5\mathcal N(4,1)$ , предложение имеет вид $y=x+\sigma\varepsilon$ с $\varepsilon\sim\mathcal N(0,1)$ . Для каждого $\sigma$ из сетки $\{0{,}2;0{,}5;1;2;4\}$ запустите восемь цепей длиной 40 000 после warm-up 10 000: старты первых четырёх равны $-4$ , остальных $+4$ , PCG64 seeds 6740–6747. Переходом между модами считайте смену знака после посещения области $|x|\ge2$ . Найдите $\sigma$ , при котором не менее четырёх цепей имеют ноль переходов, но объединённые доли $x<0$ и $x>0$ обе лежат в $[0{,}45;0{,}55]$ . Покажите восемь трасс, split- $\widehat R$ , числа переходов и объясните, почему объединённая гистограмма скрывает несмешивание.

7баллов

Используйте UCI Breast Cancer Wisconsin (Diagnostic), DOI 10.24432/C5DW2B, снимок файла wdbc.data от 1995 года: цель Diagnosis кодируйте как $y=1$ для M и $0$ для B, признаки ограничьте ровно тремя столбцами radius1, texture1, perimeter1. Отсортируйте 569 строк по ID: первые 398 образуют train, следующие 85 validation, последние 86 test; строку с пропуском удаляйте, среднее и стандартное отклонение считайте только на train. Модель: $\Pr(y_i=1)=\operatorname{logit}^{-1}(\beta_0+x_i^\mathsf T\beta)$ , все четыре коэффициента независимо имеют prior $\mathcal N(0,2^2)$ . Random-walk MH предлагает $\theta'=\theta+0{,}05z$ , $z\sim\mathcal N(0,I_4)$ ; четыре старта равны $(0,0,0,0)$ , $(2,2,2,2)$ , $(-2,-2,-2,-2)$ и $(0,2,-2,1)$ , PCG64 seeds 6751–6754, warm-up 20 000 и ещё 80 000 итераций. Нужны 95%-интервалы, split- $\widehat R<1{,}05$ , $N_{\mathrm{eff}}>400$ для каждого коэффициента и posterior predictive доля положительных на test; причинных выводов делать не следует.