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

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

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

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

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

Числитель вычислить обычно можно: это правдоподобие, умноженное на prior, — всё то, что мы уже собирали в уроке 47. А вот многомерный интеграл в знаменателе не берётся ни аналитически, ни численно: сетка из десяти узлов по каждой из двадцати координат — это $10^{20}$ точек, больше, чем секунд с начала Вселенной. Обозначим ненормированную плотность

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

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

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

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

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

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

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

Оно говорит, что в равновесии поток вероятности из $x$ в $y$ в точности компенсируется обратным потоком. Проинтегрируем по всем $x$:

$$\int\pi(dx)K(x,dy)=\pi(dy)\int K(y,dx)=\pi(dy),$$

потому что полная вероятность перейти хоть куда-нибудь равна единице. Значит, один такт цепи сохраняет целевое распределение: если $\theta_t\sim\pi$, то и $\theta_{t+1}\sim\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$. Для симметричного предложения $q(y\mid x)=q(x\mid y)$ остаётся исходная формула Метрополиса:

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

В программе всё считают в логарифмах: предложение принимается, если

$$\log u<\log\widetilde\pi(y)-\log\widetilde\pi(x), \qquad u\sim\mathrm{U}(0,1),$$

и это единственная строка, где встречается случайность помимо самого предложения.

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

Это одна фраза, и в ней весь метод. Прямое взвешивание безнадёжно: почти все взятые наугад конфигурации имеют ничтожный вес, и миллион точек несёт информации как одна. Перенос трудности из веса в способ выбора и есть изобретение 1953 года.

Хастингс сделал шаг, ради которого метод и попал в статистику: снял требование симметрии предложения. Множитель $q(x\mid y)/q(y\mid x)$ компенсирует несимметричность, и предложение можно проектировать свободно — хоть тянуть его в сторону области высокой плотности.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Отказ — это тоже такт

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

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

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

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

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

Измеряется всё это автокорреляцией функции $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}}}.$$

Смысл $\tau_{\mathrm{int}}$ прямой: дисперсия среднего по зависимой цепи равна

$$\operatorname{Var}(\bar f)\approx \frac{\tau_{\mathrm{int}}\operatorname{Var}(f)}{N} =\frac{\operatorname{Var}(f)}{N_{\mathrm{eff}}},$$

то есть цепь длины $N$ стоит ровно $N/\tau_{\mathrm{int}}$ независимых наблюдений. Сравнивать методы по длине выходного файла бессмысленно.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

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

Начните с двух мод и $\sigma=1{,}2$: посмотрите, как быстро растёт число тактов и как медленно — число переходов между модами. Затем поставьте $\sigma=0{,}05$ и убедитесь, что доля принятия почти $100\%$, а $N_{\mathrm{eff}}$ остаётся двузначным. Затем $\sigma=8$: принимается мало, трасса складывается из ступенек. Наконец, переключите цель на «узкая + широкая» — настройка, удобная для широкого холма, проскакивает узкий, а настройка для узкого еле ползёт по широкому. Одной универсальной длины шага не существует.

Отдельно поставьте старт $x=8$ при маленьком $\sigma$: первые сотни тактов цепь просто идёт к области высокой плотности, и включать их в среднее нельзя.

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

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

Диагностика $\widehat R$ сопоставляет разброс между цепями и внутри них. Если $B$ — межцепная дисперсия средних, $W$ — средняя внутрицепная дисперсия, а $n$ — длина половины цепи, то

$$\widehat{\operatorname{Var}}=\frac{n-1}{n}W+\frac{B}{n}, \qquad \widehat R=\sqrt{\frac{\widehat{\operatorname{Var}}}{W}} .$$

Эффективный размер выборки при нескольких цепях складывается: если запущено $M$ цепей длины $N$ с общим $\tau_{\mathrm{int}}$, то

$$N_{\mathrm{eff}}^{\text{всего}}=\frac{MN}{\tau_{\mathrm{int}}},$$

и четыре цепи по 25 000 тактов дают ровно столько же информации, сколько одна длиной 100 000, — зато позволяют вычислить $\widehat R$.

Значение около единицы — необходимый, но не достаточный признак. Цепи могут дружно застрять в одной моде и показать $\widehat R\approx1{,}00$, честно согласившись друг с другом в общем заблуждении.

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

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Проверка на задаче с известным ответом

Перед сложным posterior полезно устроить калибровочный стенд: взять реальные данные, но такую модель, где ответ известен точно. Возьмём SMS Spam Collection: 5574 сообщения, из них 747 помечены как спам. Модель — та же, что в уроке 47: доля спама $p$ с prior $\mathrm{Beta}(1,1)$ и биномиальным правдоподобием. Тогда posterior известен в замкнутой форме:

$$p\mid D\sim\mathrm{Beta}(1+747,\;1+4827), \qquad \mathbb E[p\mid D]=\frac{748}{5576}=0{,}13415 .$$

Теперь сделаем вид, что формулы мы не знаем, и запустим random-walk Metropolis по $p$ с $\sigma=0{,}02$: 120 000 тактов, первые 20 000 в прогрев. Цепь приняла $27{,}2\%$ предложений, дала среднее $0{,}13410$, $\tau_{\mathrm{int}}=5{,}8$ и $N_{\mathrm{eff}}=8645$.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

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

Геометрия posterior: реальные данные

Возьмём набор Breast Cancer Wisconsin из sklearn: 569 наблюдений, из них 212 злокачественных. Оставим два признака — mean radius и mean perimeter — и заметим главное: у почти круглого пятна периметр пропорционален радиусу, поэтому корреляция признаков равна $0{,}998$. Модель:

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

а ненормированный логарифм posterior равен

$$\log\widetilde\pi(\beta) =\sum_i\Bigl[y_i z_i-\log\bigl(1+e^{z_i}\bigr)\Bigr] -\frac{\|\beta\|^2}{2\cdot4}, \qquad z_i=\beta_0+x_i^\top\beta .$$

Почти дублирующие признаки делают posterior вытянутым: корреляция коэффициентов равна $-0{,}963$, отношение осей эллипса $11{,}9$. Изотропное предложение обязано быть узким по короткой оси — и потому ползёт вдоль длинной.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Интересно, что́ именно posterior знает уверенно. Отдельные коэффициенты почти не определены: $\beta_{\text{radius}}=-2{,}89$ с интервалом $[-5{,}26;\,-0{,}50]$ шириной $4{,}76$, а $\beta_{\text{perimeter}}=6{,}69$ с интервалом $[4{,}26;\,9{,}10]$. Зато их сумма определена почти вчетверо точнее: $3{,}80$ с интервалом $[3{,}17;\,4{,}47]$ шириной $1{,}30$.

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

Плотность часто записывают в физической форме

$$\pi(x)\propto e^{-E(x)/T}, \qquad E(x)=-T\log\widetilde\pi(x).$$

Энергия $E$ мала в вероятных состояниях, температура $T$ сглаживает рельеф. При $T\to0$ пики сжимаются к точкам максимума, при больших $T$ рельеф становится почти плоским. Отношение принятия принимает вид

$$\frac{\widetilde\pi(y)}{\widetilde\pi(x)}=e^{-\Delta E/T}, \qquad \Delta E=E(y)-E(x),$$

а вероятность за один шаг пересечь барьер высотой $\Delta E$ падает экспоненциально. Отсюда и берутся нулевые переходы на рис. 67.3.

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

Сколько стоит одна честная цифра

Оценив среднее, надо честно сообщить его погрешность. Наивная формула из урока 46 предполагает независимость и потому в MCMC врёт. Правильная версия — Monte Carlo standard error:

$$\operatorname{MCSE}(\bar f) \approx\sqrt{\frac{\widehat{\operatorname{Var}}(f)}{N_{\mathrm{eff}}}} =\sqrt{\frac{\tau_{\mathrm{int}}\,\widehat{\operatorname{Var}}(f)}{N}} .$$

Проверить её можно так же, как мы проверяли интервалы в уроке 46: многократно повторить всю процедуру и посчитать, как часто интервал накрывает известное истинное значение.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Отсюда простое правило работы: любое число из цепи публикуется вместе с $N_{\mathrm{eff}}$ и MCSE, а количество значащих цифр в ответе определяется MCSE, а не тем, сколько их напечатал компьютер.

Другая ветвь Монте-Карло: Соболь и равномерность вместо случайности

У задачи «оценить интеграл» есть и второй путь, и его в 1967 году проложил Илья Меерович Соболь. В работе «О распределении точек в кубе и приближённом вычислении интегралов» он построил ЛП$\tau$-последовательности — не случайные, а нарочито равномерные наборы точек, у которых каждый двоичный блок координат заполняется без сгущений и пустот. Для гладких функций ошибка квази-Монте-Карло убывает почти как $1/N$ вместо $1/\sqrt N$:

$$\bigl|\widehat I_N-I\bigr|\le V(f)\,D^{*}_N, \qquad D^{*}_N=O\!\left(\frac{(\log N)^{d}}{N}\right),$$

где $V(f)$ — вариация функции, а $D^{*}_N$ — дискрепанс, то есть отклонение набора точек от идеальной равномерности. Это неравенство Коксмы–Хлавки: детерминированная оценка, без всяких «с вероятностью $0{,}95$».

Две ветви решают разные задачи, и это стоит помнить. ЛП$\tau$ великолепен для гладкого интеграла умеренной размерности по известной области; MCMC незаменим там, где область — это сама неизвестная масса вероятности в двадцатимерном пространстве, и найти её можно, только ползая по отношениям плотностей. Советская школа вычислительной вероятности разрабатывала обе ветви: Соболь и Ермаков строили теорию Монте-Карло для интегралов и интегральных уравнений, тогда как физики пользовались метрополисовскими цепями. Сегодняшние библиотеки содержат обе кнопки, и выбирать между ними — часть работы.

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

Успешная MCMC-работа состоит из трёх слоёв.

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

Второй — вычислительная диагностика: несколько цепей из разнесённых стартов, трассы, split-$\widehat R$, $\tau_{\mathrm{int}}$, $N_{\mathrm{eff}}$, MCSE. Она отвечает на вопрос, успела ли цепь исследовать существенные области за отведённое время. Доля принятия относится сюда — и только сюда.

Третий — предметная проверка: осмысленна ли сама вероятностная модель, откуда взяты данные, что означает коэффициент. Идеальная цепь на плохой модели даёт уверенно неверный ответ, и никакой $\widehat R$ этого не заметит.

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

Задачи