Введение

Алгоритм Монте-Карло

В статистике, выборка Гиббса или сэмплинг Гиббса – это алгоритм цепи Маркова Монте-Карло (MCMC) для выборки из заданного многомерного распределения вероятностей, когда непосредственная выборка из совместного распределения затруднена, но выборка из условного распределения более целесообразна. Эта последовательность может быть использована для аппроксимации совместного распределения (например, для построения гистограммы распределения); для аппроксимации краевого распределения одной из переменных или некоторого подмножества переменных (например, неизвестных параметров или латентных переменных); или для вычисления интеграла (например, математического ожидания одной из переменных). Как правило, некоторые переменные соответствуют наблюдениям, значения которых известны и, следовательно, не требуют выборки. Выборка Гиббса часто используется как средство статистического вывода, особенно байесовского вывода. Это рандомизированный алгоритм (то есть алгоритм, использующий случайные числа) и является альтернативой детерминированным алгоритмам для статистического вывода, таким как алгоритм максимизации ожиданий (EM). Как и другие алгоритмы MCMC, выборка Гиббса генерирует цепь Маркова выборок, каждая из которых коррелирует с соседними выборками. В результате, необходимо соблюдать осторожность, если требуются независимые выборки. Обычно выборки с начала цепи (период прогрева) могут неточно представлять желаемое распределение и обычно отбрасываются.

Введение

Пробоиночный метод Гиббса назван в честь физика Джозии Вилларда Гиббса, в связи с аналогией между алгоритмом выборки и статистической физикой. Алгоритм был описан братьями Стюартом и Дональдом Геманами в 1984 году, примерно через восемь десятилетий после смерти Гиббса, и получил широкое распространение в статистическом сообществе для вычисления предельных распределений вероятностей, особенно апостериорного распределения. В своей базовой версии выборка Гиббса является частным случаем алгоритма Метрополиса — Хэстингса. Однако в расширенных версиях (см. ниже) его можно рассматривать как общую структуру для выборки из большого набора переменных путем последовательной выборки каждой переменной (или, в некоторых случаях, каждой группы переменных). Он может включать алгоритм Метрополиса — Хэстингса (или методы, такие как метод нарезки) для реализации одного или нескольких шагов выборки. Выборка Гиббса применима, когда совместное распределение неизвестно явно или сложно для непосредственной выборки, но условное распределение каждой переменной известно и легко (или, по крайней мере, проще) для выборки. Алгоритм выборки Гиббса генерирует выборку из распределения каждой переменной последовательно, при условии текущих значений остальных переменных. Можно показать, что последовательность выборок образует цепь Маркова, а стационарное распределение этой цепи Маркова является искомым совместным распределением. Выборка Гиббса особенно хорошо подходит для выборки апостериорного распределения в байесовской сети, поскольку байесовские сети обычно задаются как набор условных распределений.

Реализация

Выборка Гиббса, в своей базовой реализации, является частным случаем алгоритма Метрополиса–Хестингса. Преимущество выборки Гиббса заключается в том, что для многомерного распределения проще проводить выборку из условных распределений, чем вычислять маргинальное распределение путем интегрирования по совместному распределению. Предположим, мы хотим получить *n* образцов из совместного распределения. Обозначим *i*-й образец как Мы действуем следующим образом:

Начинаем с некоторого начального значения. Нам нужен следующий образец. Обозначим этот образец как . Поскольку является вектором, мы выбираем каждый компонент этого вектора, , из распределения этого компонента, обусловленного всеми остальными компонентами, уже выбранными на данный момент. Однако есть нюанс: мы обуславливаемся компонентами до , а затем обуславливаемся компонентами , начиная с . Чтобы этого добиться, мы выбираем компоненты последовательно, начиная с первого компонента. Более формально, для выборки , мы обновляем его в соответствии с распределением, заданным . Мы используем значение, которое *i*-й компонент имел в (*i*-1)-м образце, а не в *i*-м образце. Повторяем описанные выше шаги *n* раз.

Вывод

Выборка Гиббса обычно используется для статистического вывода (например, для определения наилучшего значения параметра, такого как определение количества людей, которые, вероятно, совершат покупку в определенном магазине в определенный день, кандидата, за которого избиратель, скорее всего, проголосует, и т. д.). Идея заключается в том, чтобы включить наблюдаемые данные в процесс выборки, создавая отдельные переменные для каждого фрагмента наблюдаемых данных и фиксируя соответствующие переменные на их наблюдаемых значениях, а не выполняя выборку из этих переменных. Распределение оставшихся переменных, таким образом, фактически является апостериорным распределением, обусловленным наблюдаемыми данными. Наиболее вероятное значение желаемого параметра (мода) можно просто выбрать, выбрав значение выборки, которое встречается наиболее часто; это по существу эквивалентно оценке апостериорной максимальной правдоподобности параметра. (Поскольку параметры обычно непрерывны, часто необходимо "разбить" отобранные значения на один из конечного числа диапазонов или "интервалов", чтобы получить значимую оценку моды.) Однако чаще всего выбирается ожидаемое значение (среднее или среднее арифметическое) отобранных значений; это байесовский оценщик, который использует дополнительные данные обо всем распределении, доступные из байесовской выборки, в то время как алгоритм максимизации, такой как максимизация ожидания (EM), способен возвращать только одну точку из распределения. Например, для унимодального распределения среднее значение (ожидаемое значение) обычно близко к моде (наиболее частому значению), но если распределение скошено в одном направлении, среднее значение будет смещено в этом направлении, что эффективно учитывает дополнительную массу вероятности в этом направлении. (Если распределение мультимодальное, ожидаемое значение может не возвращать значимую точку, и любая из мод, как правило, является лучшим выбором.) Хотя некоторые из переменных обычно соответствуют интересующим параметрам, другие являются неинтересными ("побочными") переменными, введенными в модель для правильного выражения взаимосвязей между переменными. Хотя отобранные значения представляют собой совместное распределение по всем переменным, побочные переменные можно просто игнорировать при вычислении ожидаемых значений или мод; это эквивалентно маргинализации по побочным переменным. Когда требуется значение для нескольких переменных, ожидаемое значение просто вычисляется для каждой переменной отдельно. (Однако при вычислении моды все переменные должны рассматриваться вместе.) Обучение с учителем, обучение без учителя и полу-контролируемое обучение (также известное как обучение с пропущенными значениями) можно обработать, просто зафиксировав значения всех переменных, значения которых известны, и выполняя выборку из остальных. Для наблюдаемых данных будет одна переменная для каждого наблюдения, а не, например, одна переменная, соответствующая среднему значению выборки или дисперсии выборки набора наблюдений. Фактически, обычно не будет переменных, соответствующих таким понятиям, как "среднее значение выборки" или "дисперсия выборки". Вместо этого в таком случае будут переменные, представляющие неизвестное истинное среднее и истинную дисперсию, и определение значений выборки для этих переменных автоматически происходит в результате работы сэмплера Гиббса. Обобщенные линейные модели (т. е. варианты линейной регрессии) иногда также можно обработать с помощью выборки Гиббса. Например, пробит-регрессия для определения вероятности данного бинарного (да/нет) выбора, с нормальными априорными распределениями, заданными для коэффициентов регрессии, может быть реализована с помощью выборки Гиббса, поскольку можно добавить дополнительные переменные и воспользоваться свойством сопряженности. Однако логистическую регрессию таким образом обработать нельзя. Один из вариантов — аппроксимировать логистическую функцию смесью (обычно 7–9) нормальных распределений. Однако чаще всего вместо выборки Гиббса используется метод Метрополиса — Хэстингса.

Колапс Дирихлетовских распределений

В иерархических байесовских моделях с категориальными переменными, таких как скрытое распределение Дирихле и различные другие модели, используемые в обработке естественного языка, часто применяется исключение (коллапс) распределений Дирихле, которые обычно используются в качестве априорных распределений для категориальных переменных. Результатом этого исключения становится введение зависимостей между всеми категориальными переменными, зависящими от заданного распределения Дирихле, а совместное распределение этих переменных после исключения представляет собой распределение Дирихле-мультиномиальное. Условное распределение данной категориальной переменной в этом распределении, при условии остальных, принимает чрезвычайно простую форму, что облегчает выборку Гиббса по сравнению со случаем, когда исключение не выполнялось. Правила таковы: исключение априорного узла Дирихле влияет только на родительский и дочерние узлы этого априорного распределения. Поскольку родительский узел часто является константой, обычно беспокоиться нужно только о дочерних узлах. Исключение априорного распределения Дирихле вводит зависимости между всеми категориальными дочерними узлами, зависящими от этого априорного распределения, но не вводит дополнительных зависимостей между другими категориальными дочерними узлами. (Это важно учитывать, например, когда существует несколько априорных распределений Дирихле, связанных одним и тем же гиперприором. Каждое априорное распределение Дирихле можно исключать независимо, и оно влияет только на его прямых дочерних узлов.) После исключения условное распределение одного зависимого дочернего узла при условии остальных принимает очень простую форму: вероятность наблюдения данного значения пропорциональна сумме соответствующего гиперприора для этого значения и количеству всех остальных зависимых узлов, принимающих то же значение. Узлы, не зависящие от одного и того же априорного распределения, не должны учитываться. То же правило применяется и в других итеративных методах вывода, таких как вариационный байесовский метод или максимизация ожидания; однако, если метод предполагает сохранение частичных подсчетов, то частичный подсчет для рассматриваемого значения должен быть суммирован по всем остальным зависимым узлам. Иногда эта суммарная частичная величина называется ожидаемым подсчетом или аналогичным образом. Вероятность пропорциональна полученному значению; фактическая вероятность должна быть определена путем нормализации относительно всех возможных значений, которые может принимать категориальная переменная (то есть, путем суммирования вычисленного результата для каждого возможного значения категориальной переменной и деления всех вычисленных результатов на эту сумму). Если данный категориальный узел имеет зависимые дочерние узлы (например, когда это скрытая переменная в смесительной модели), значение, вычисленное на предыдущем шаге (ожидаемый подсчет плюс априорное распределение, или что бы то ни было вычисленное), должно быть умножено на фактические условные вероятности (а не на вычисленное значение, пропорциональное вероятности!) всех дочерних узлов, обусловленные их родительскими узлами. См. статью о распределении Дирихле-мультиномиальное для подробного обсуждения. В случае, когда групповая принадлежность узлов, зависящих от заданного априорного распределения Дирихле, может динамически изменяться в зависимости от другой переменной (например, категориальной переменной, индексированной другой скрытой категориальной переменной, как в тематической модели), те же ожидаемые подсчеты все еще вычисляются, но необходимо делать это осторожно, чтобы был включен правильный набор переменных. См. статью о распределении Дирихле-мультиномиальное для более подробного обсуждения, в том числе в контексте тематической модели.

Срыв других приоров

В общем, любой сопряженный априорный закон можно исключить, если у его единственных потомков распределения сопряжены с ним. Соответствующая математика обсуждается в статье о составных распределениях. Если есть только один дочерний узел, результат часто предполагает известное распределение. Например, исключение дисперсии, распределенной по обратному гамма-закону, из сети с одним гауссовским потомком даст t-распределение Стьюдента. (Более того, исключение среднего и дисперсии одного гауссовского потомка также даст t-распределение Стьюдента, при условии, что оба сопряжены, то есть гауссовское среднее, дисперсия, распределенная по обратному гамма-закону.) Если есть несколько дочерних узлов, они все станут зависимыми, как в случае Дирихле-категориального распределения. Полученное совместное распределение будет иметь замкнутую форму, в некотором смысле напоминающую составное распределение, хотя оно будет представлять собой произведение ряда факторов, по одному для каждого дочернего узла. Кроме того, и что наиболее важно, полученное условное распределение одного из дочерних узлов при заданных остальных (а также при заданных родителях исключаемого узла, но не при заданных потомках дочерних узлов) будет иметь ту же плотность, что и апостериорное предсказательное распределение всех оставшихся дочерних узлов. Более того, апостериорное предсказательное распределение имеет ту же плотность, что и базовое составное распределение одного узла, хотя и с другими параметрами. Общая формула приведена в статье о составных распределениях. Например, для байесовской сети с набором условно независимых, одинаково распределенных гауссовских узлов с сопряженными априорными законами, заданными для среднего и дисперсии, условное распределение одного узла при заданных остальных после исключения как среднего, так и дисперсии будет t-распределением Стьюдента. Аналогично, результат исключения гамма-априорного закона для ряда узлов, распределенных по закону Пуассона, приводит к тому, что условное распределение одного узла при заданных остальных принимает вид отрицательного биномиального распределения. В этих случаях, когда исключение приводит к хорошо известному распределению, часто существуют эффективные процедуры выборки, и их использование часто (хотя и не обязательно) будет более эффективным, чем не исключать, а вместо этого выбирать как априорные, так и дочерние узлы по отдельности. Однако в случае, когда составное распределение не является хорошо известным, из него может быть нелегко проводить выборку, поскольку оно обычно не будет принадлежать к экспоненциальному семейству и, как правило, не будет лог-вогнутым (что облегчило бы выборку с использованием адаптивной выборки с отбраковкой, поскольку замкнутая форма всегда существует). В случае, когда дочерние узлы исключаемых узлов сами имеют потомков, условное распределение одного из этих дочерних узлов при заданных всех остальных узлах в графе должно учитывать распределение этих потомков второго уровня. В частности, полученное условное распределение будет пропорционально произведению составного распределения, как определено выше, и условных распределений всех дочерних узлов при заданных их родителях (но не при заданных их собственных потомках). Это следует из того факта, что полное условное распределение пропорционально совместному распределению. Если дочерние узлы исключаемых узлов непрерывны, это распределение, как правило, не будет иметь известной формы и может быть трудно для выборки, несмотря на то, что замкнутую форму можно записать, по тем же причинам, что и описанные выше для плохо известных составных распределений. Однако в том случае, если дочерние узлы дискретны, выборка возможна, независимо от того, являются ли потомки этих дочерних узлов непрерывными или дискретными. Фактически, принцип, задействованный здесь, достаточно подробно описан в статье о распределении Дирихле-многочлена.

Пробоотборник Гиббса с заказным перерасслаблением

Пробоотборник Гиббса с упорядоченной перерасслаблением на каждом шаге генерирует заданное нечетное число кандидатов для и сортирует их вместе с единственным значением для в соответствии с некоторым четко определенным порядком. Если является -м наименьшим в отсортированном списке, то выбирается как -й наибольший в отсортированном списке. Более подробную информацию можно найти в работе Neal (1995).

Другие расширения

Также возможно расширить выборку Гиббса различными способами. Например, в случае переменных, из которых сложно проводить выборку, учитывая их условное распределение, можно использовать одну итерацию алгоритма выборок на срезах (slice sampling) или алгоритма Метрополиса — Гастингса для выборки значений этих переменных. Также можно включать переменные, которые не являются случайными, но чьи значения детерминированно вычисляются на основе других переменных. Обобщенные линейные модели, например, логистическая регрессия (также известные как "модели максимальной энтропии"), могут быть включены таким образом. (Например, BUGS допускает подобное комбинирование моделей.)

Режимы отказов

Есть два способа, которыми выборка Гиббса может завершиться неудачей. Первый – это когда существуют изолированные области состояний с высокой вероятностью, между которыми нет переходов. Например, рассмотрим распределение вероятностей по 2-битным векторам, где векторы (0,0) и (1,1) каждый имеют вероятность ½, а два других вектора (0,1) и (1,0) имеют вероятность 0. Выборка Гиббса застрянет в одном из двух векторов с высокой вероятностью и никогда не достигнет другого. В более общем случае, для любого распределения по многомерным векторам с вещественными значениями, если два конкретных элемента вектора идеально коррелированы (или идеально антикоррелированы), эти два элемента зафиксируются, и выборка Гиббса не сможет их изменить. Вторая проблема может возникнуть даже тогда, когда все состояния имеют ненулевую вероятность и существует только одна изолированная область состояний с высокой вероятностью. Например, рассмотрим распределение вероятностей по 100-битным векторам, где вектор, состоящий только из нулей, встречается с вероятностью ½, а все остальные векторы равновероятны, то есть каждый имеет вероятность. Если вы хотите оценить вероятность нулевого вектора, было бы достаточно взять 100 или 1000 выборок из истинного распределения. Это, скорее всего, даст результат, очень близкий к ½. Но для получения того же результата вам, вероятно, потребуется взять более чем выборок из выборки Гиббса. Ни один компьютер не сможет выполнить это за всю жизнь. Эта проблема возникает независимо от длительности периода прогрева. Это связано с тем, что в истинном распределении нулевой вектор встречается в половине случаев, и эти случаи случайным образом перемешаны с ненулевыми векторами. Даже небольшая выборка покажет как нулевые, так и ненулевые векторы. Однако выборка Гиббса будет чередоваться, возвращая только нулевой вектор в течение длительных периодов (примерно в ряд), а затем только ненулевые векторы в течение длительных периодов (примерно в ряд). Таким образом, сходимость к истинному распределению происходит крайне медленно, требуя гораздо больше, чем шагов; выполнение такого количества шагов не представляется вычислительно возможным за разумное время. Медленная сходимость здесь может рассматриваться как следствие проклятия размерности. Эту проблему можно решить, выполнив блочную выборку всего 100-битного вектора сразу. (Это предполагает, что 100-битный вектор является частью большего набора переменных. Если этот вектор является единственным объектом выборки, то блочная выборка эквивалентна полному отсутствию выборки Гиббса, что, по гипотезе, было бы затруднительно.)

Программное обеспечение

Программное обеспечение OpenBUGS (Bayesian inference Using Gibbs Sampling) выполняет байесовский анализ сложных статистических моделей, используя метод Монте-Карло на цепях Маркова. JAGS (Just another Gibbs sampler) – это программа, распространяемая по лицензии GPL, для анализа байесовских иерархических моделей с использованием метода Монте-Карло на цепях Маркова. Church – это свободное программное обеспечение для выполнения вывода Гиббса для произвольных распределений, заданных в виде вероятностных программ. PyMC – это библиотека Python с открытым исходным кодом для байесовского обучения обобщенным вероятностным графическим моделям. Turing – это библиотека Julia с открытым исходным кодом для байесовского вывода, использующая вероятностное программирование.