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

03.10.2022

Семплирование по Гиббсу — алгоритм для генерации выборки совместного распределения множества случайных величин. Он используется для оценки совместного распределения и для вычисления интегралов методом Монте-Карло. Этот алгоритм является частным случаем алгоритма Метрополиса-Гастингса и назван в честь физика Джозайи Гиббса.

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

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

Алгоритм

Пусть есть совместное распределение p ( x 1 , . . . , x d ) {displaystyle p(x_{1},...,x_{d})} для d {displaystyle d} случайных величин, причём d {displaystyle d} может быть очень большим. Пусть на шаге t {displaystyle t} мы уже выбрали какое-то значение X = { x i t } {displaystyle X={x_{i}^{t}}} . На каждом шаге делаются следующие действия:

  • Выбирается индекс i : ( 1 ≤ i ≤ d {displaystyle i:(1leq ileq d} ).
  • x i t + 1 {displaystyle x_{i}^{t+1}} выбирается по распределению p ( x i | x 1 t , . . . , x i − 1 t , x i + 1 t , . . . , x d t ) {displaystyle p(x_{i}|x_{1}^{t},...,x_{i-1}^{t},x_{i+1}^{t},...,x_{d}^{t})} , а для остальных индексов значение не меняется: x j t + 1 = x j t {displaystyle x_{j}^{t+1}=x_{j}^{t}} (j≠i).

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

    Пример

    Пусть есть совместное распределение p ( x 1 , x 2 , x 3 ) {displaystyle p(x_{1},x_{2},x_{3})} из трех случайных величин, каждая из которых находится в диапазоне от 0 до 10.

    Примем, что первоначальное значение вектора, от которого начнется итерационный процесс, будет X = { 5 , 2 , 7 } {displaystyle X={5,2,7}} .

    Далее фиксируем x 2 {displaystyle x_{2}} и x 3 {displaystyle x_{3}} , после чего рассчитываем по известной заранее формуле условную вероятность p ( x 1 | x 2 , x 3 ) {displaystyle p(x_{1}|x_{2},x_{3})} , то есть p ( x 1 | x 2 = 2 , x 3 = 7 ) {displaystyle p(x_{1}|x_{2}=2,x_{3}=7)} , получая некоторый график плотности вероятности от переменной x 1 {displaystyle x_{1}} . То, что изначально x 1 {displaystyle x_{1}} мы положили равным 5, забываем, больше это значение не понадобится.

    Теперь необходимо выполнить семплирование — сгенерировать новое случайное значение для x 1 {displaystyle x_{1}} в соответствии с полученной плотностью вероятности. Семплирование можно сделать, например, по алгоритму выборки с отклонением. Для этого генерируется случайное число с равномерным распределением от 0 до 10, после чего для этого сгенерированного числа вычисляется его вероятность по графику плотности вероятности p ( x 1 | x 2 = 2 , x 3 = 7 ) {displaystyle p(x_{1}|x_{2}=2,x_{3}=7)} .

    Например, пусть сгенерировалось случайное число 4 и по графику плотности его вероятность равна 0.2. Тогда, в соответствии с алгоритмом выборки с отклонением, мы принимаем это сгенерированное число с вероятностью 0.2. А для этого, в свою очередь, генерируем ещё одно случайное число от 0 до 1 с равномерным распределением, и, если сгенерировалось число меньше 0.2, то мы принимаем число 4 как успешное. Иначе повторяем сначала — генерируем ещё одно число (например выпадает 3), для него находим вероятность (например, 0.3), для него генерируем ещё число от 0 до 1 (например, 0.1) и тогда уже принимаем окончательно, что на этой итерации x 1 = 3 {displaystyle x_{1}=3} .

    Далее необходимо повторить все действия выше с величиной x 2 {displaystyle x_{2}} , причём x 1 {displaystyle x_{1}} мы уже используем «новое» — в нашем примере равное 3. Так, рассчитываем плотность вероятности p ( x 2 | x 1 = 3 , x 3 = 7 ) {displaystyle p(x_{2}|x_{1}=3,x_{3}=7)} , генерируем снова случайное число на роль кандидата нового значения x 2 {displaystyle x_{2}} , делаем выборку с отклонением и повторяем её в случае, если значение «отклонено».

    Аналогично действия повторяются для x 3 {displaystyle x_{3}} с новыми значениями x 1 {displaystyle x_{1}} и x 2 {displaystyle x_{2}} . Первая итерация алгоритма семплирования по Гиббсу завершена. Через несколько сотен/тысяч таких итераций случайные значения должны прийти к максимуму своей плотности, который может быть расположен достаточно далеко от нашего первого приближения X = { 5 , 2 , 7 } {displaystyle X={5,2,7}} и семплироваться в той области. Дальнейшая тысяча итераций может уже использоваться по назначению (для поиска математического ожидания, например) как образец значений искомого распределения, не зависящих от первоначального вектора X = { 5 , 2 , 7 } {displaystyle X={5,2,7}} .