Сэмплирование Гиббса

Сэмплирование Гиббса — это алгоритм марковского цепного Монте-Карло (MCMC) для получения последовательности наблюдений, аппроксимирующих многомерное распределение вероятностей. Он применяется, когда прямое семплирование из совместного распределения затруднено, но возможно семплирование из каждого условного распределения переменной при фиксированных значениях остальных. Метод назван в честь американского физика Джозайи Уилларда Гиббса, хотя его современная формулировка в контексте байесовской статистики и вычислительной физики была разработана в 1980-х годах.

История

Предпосылки и ранние работы

Методы MCMC, включая сэмплирование Гиббса, берут начало в работах Станислава Улама, Джона фон Неймана и Николаса Метрополиса в 1940-х — 1950-х годах в Лос-Аламосской национальной лаборатории. Первоначально они использовались для моделирования нейтронных потоков в ядерных реакторах. Алгоритм Метрополиса — Гастингса, опубликованный в 1953 году, стал основой для многих последующих методов.

Формализация Гиббс-сэмплера

В 1984 году Стюарт Герман и Дональд Джеман в статье «Стохастическая релаксация, распределения Гиббса и байесовская реставрация изображений» впервые применили итеративное семплирование из условных распределений для восстановления изображений. Они назвали этот процесс «сэмплированием Гиббса» из-за его связи с распределением Гиббса из статистической механики. Позднее, в 1990-х годах, метод был популяризирован в байесовской статистике благодаря работам Адриана Смита, Алана Гельмана и других.

Алгоритм

Основная идея

Пусть требуется получить выборку из многомерного распределения \( P(x_1, x_2, \dots, x_n) \). Сэмплирование Гиббса строит цепь Маркова, последовательно обновляя каждую переменную. На каждом шаге \( t \) для каждой переменной \( x_i \) новое значение \( x_i^{(t+1)} \) семплируется из условного распределения \( P(x_i \mid x_1^{(t+1)}, \dots, x_{i-1}^{(t+1)}, x_{i+1}^{(t)}, \dots, x_n^{(t)}) \), то есть при фиксированных текущих значениях всех остальных переменных.

Пошаговое описание

1. Инициализация: выбираются начальные значения \( x_1^{(0)}, x_2^{(0)}, \dots, x_n^{(0)} \).
2. Итерация для \( t = 0, 1, 2, \dots \):

• Семплировать \( x_1^{(t+1)} \sim P(x_1 \mid x_2^{(t)}, x_3^{(t)}, \dots, x_n^{(t)}) \).
• Семплировать \( x_2^{(t+1)} \sim P(x_2 \mid x_1^{(t+1)}, x_3^{(t)}, \dots, x_n^{(t)}) \).
• ...
• Семплировать \( x_n^{(t+1)} \sim P(x_n \mid x_1^{(t+1)}, x_2^{(t+1)}, \dots, x_{n-1}^{(t+1)}) \).

1. Повторение до сходимости цепи.

Особенности

• Сходимость: при выполнении условий эргодичности (неприводимость, апериодичность) цепь Маркова сходится к целевому распределению. Первые \( B \) итераций (burn-in) обычно отбрасываются для уменьшения влияния начальных значений.
• Корреляция: последовательные сэмплы автокоррелированы, поэтому для получения независимых наблюдений может потребоваться прореживание (thinning) — сохранение каждого \( k \)-го сэмпла.
• Стохастическая релаксация: в оригинальной работе Германа и Джемана алгоритм рассматривался как метод стохастической релаксации, где переменные обновляются случайным образом, а не детерминированно.

Применение

Байесовская статистика

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

• Латентное распределение Дирихле (LDA): тематическое моделирование текстов, где сэмплирование Гиббса применяется для оценки распределений тем по словам и документам.
• Модели линейной регрессии: оценка параметров регрессии и дисперсии при априорных распределениях.
• Гауссовские процессы: семплирование гиперпараметров в моделях гауссовских процессов.

Машинное обучение

• Марковские случайные поля (MRF): сэмплирование Гиббса используется для вывода в графовых вероятностных моделях, например, в реставрации изображений или сегментации.
• Глубокие порождающие модели: в ограниченных машинах Больцмана (RBM) и глубоких сетях доверия (DBN) сэмплирование Гиббса применяется для обучения и генерации данных.

Вычислительная физика и химия

• Моделирование спиновых систем: например, модель Изинга, где сэмплирование Гиббса используется для вычисления средних намагниченностей и корреляций.
• Молекулярная динамика: оценка конформаций молекул и свободной энергии.

Обработка изображений

• Восстановление изображений: удаление шума и восстановление пропущенных пикселей на основе модели Марковского случайного поля.
• Сегментация: классификация пикселей по классам с использованием априорных распределений.

Варианты и модификации

Блочное сэмплирование Гиббса

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

Сэмплирование Гиббса с коллапсом

В некоторых моделях часть переменных может быть проинтегрирована аналитически, что уменьшает размерность пространства. Например, в LDA часто используется коллапсированное сэмплирование Гиббса, где интегрируются параметры распределений тем.

Сэмплирование Гиббса с метрополисом-в-гиббсе

Если условное распределение не является стандартным, для семплирования из него может использоваться шаг алгоритма Метрополиса — Гастингса. Это позволяет применять метод к более широкому классу моделей.

Адаптивное сэмплирование Гиббса

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

Ограничения и критика

Высокая корреляция сэмплов

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

Проблемы сходимости

• Медленная сходимость: при сильной корреляции между переменными (например, в моделях с высокой размерностью) цепь может «застревать» в определённых областях пространства.
• Неприводимость: если условные распределения не позволяют достичь всех областей пространства, цепь может быть неэргодичной.

Чувствительность к инициализации

Начальные значения могут существенно влиять на время сходимости. Неудачный выбор может привести к долгому периоду burn-in или даже к неправильной аппроксимации.

Вычислительная сложность

На каждом шаге требуется вычисление условных распределений, что может быть затратно для моделей с большим числом переменных или сложными зависимостями. Для высокоразмерных задач (например, с тысячами переменных) сэмплирование Гиббса может быть непрактичным.

Сравнение с другими методами MCMC

Алгоритм Метрополиса — Гастингса

• Гибкость: Метрополис — Гастингс может работать с любым целевым распределением, не требуя аналитического вида условных распределений.
• Эффективность: сэмплирование Гиббса часто более эффективно, если условные распределения легко семплируются, так как оно не требует этапа принятия/отклонения.

Гамильтоново сэмплирование (HMC)

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

Сэмплирование срезов (Slice sampling)

• Преимущества: автоматически адаптирует масштаб шага, не требуя ручной настройки.
• Недостатки: может быть медленным для многомерных распределений.

Примеры

Пример 1: Двумерное нормальное распределение

Пусть \( (X, Y) \) имеют совместное нормальное распределение с нулевыми средними, единичными дисперсиями и корреляцией \( \rho \). Условные распределения:

• \( X \mid Y \sim \mathcal{N}(\rho Y, 1 - \rho^2) \)
• \( Y \mid X \sim \mathcal{N}(\rho X, 1 - \rho^2) \)

Алгоритм:

1. Инициализировать \( X^{(0)} = 0, Y^{(0)} = 0 \).
2. На каждом шаге \( t \):

• Семплировать \( X^{(t+1)} \sim \mathcal{N}(\rho Y^{(t)}, 1 - \rho^2) \).
• Семплировать \( Y^{(t+1)} \sim \mathcal{N}(\rho X^{(t+1)}, 1 - \rho^2) \).

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

Пример 2: Латентное распределение Дирихле (LDA)

В тематическом моделировании LDA сэмплирование Гиббса используется для оценки распределения тем по словам. Для каждого документа и каждого слова в нём семплируется тема из условного распределения, зависящего от текущих назначений тем в других документах. Этот процесс повторяется до сходимости, после чего оцениваются параметры модели.

Интересные факты

• Название «сэмплирование Гиббса» не связано напрямую с работами Джозайи Уилларда Гиббса по статистической механике, а было выбрано авторами метода из-за математической аналогии с распределением Гиббса.
• В 1990-х годах метод стал одним из ключевых инструментов в байесовской статистике, что привело к бурному развитию вычислительных методов в этой области.
• Сэмплирование Гиббса является частным случаем алгоритма Метрополиса — Гастингса, где предложения для каждой переменной всегда принимаются с вероятностью 1.

Источники

• Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence.
• Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.
• Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods (2nd ed.). Springer.
• Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer.
• Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press.