Методы Монте Карло по схеме Марковских цепей

Идея метода MCMC заключается в запуске алгоритма Марковской цепи, которая обладает априорным распределением p0(T). Общий алгоритм цепи заключается в последовательном совершении переходов от одного распределения плотности вероятности к другому. В результате переходов мы получаем выборку из большого числа распределений. Вероятности перехода в момент времени t можно описать, как qr (Tr+1|Tt), а выборка генерируется следующим способом:

T1 p0(T),.

T2 q1(T2|T1),.

Tr qr?1(Tr|Tr?1) P. E. Rossi, G. M. Allenby and R. McCulloch Bayesian Statistics and Marketing 2005 John Wiley & Sons, Ltd. ISBN: 0−470−86 367−6 P 50−51.

Заметим, что при таком подходе распределение на шаге t всегда зависит от распределения на шаге t-1. Следовательно, генерируемая выборка не является набором независимых случайных величин. Тем не менее, для conjoint анализа нам необходимо получить выборку из независимых величин. Для достижения независимости между величинами необходимо проредить полученную выборку T1,. .. , TT., сохраняя лишь каждый m-ый отсчет, где m достаточно велико Они же.

Работая с цепями Маркова, исследователи-маркетологи обращают должное внимание на несколько ограничений Марковской цепи.

• · Проверка эргодичности Марковской цепи;
• · Проверка сходимости и burn in Марковской цепи.

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

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

Допустим, мы хотим убедиться, что критерий эргодичности выполняется. Сперва мы должны удостовериться в том, что Марковская цепь может считаться однородной (гомогенной) если вероятность перехода qt (Tr+1|Tr) не зависит от момента времени n; вероятности переходов неизменны от шага к шагу и описываются следующим равенством:

qt (Tr+1|Tr) = q (Tr+1|Tr).

Допустим, мы будем работать дальше только с однородными Марковскими цепями. Мы создали набор значений однородной цепи при условии распределения p0(T) и вероятности перехода qt (Tr+1|Tr) = q (Tr+1|Tr). При некоторых ограничениях для функции условной вероятности и при r стремящемся к бесконечности, распределение (Tr|Tr-0) сойдется в устойчивое уникальное распределение. Такое распределение принято называть инвариантным.

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

Также мы можем записать уравнение для оценки апостериорного распределения для любой функции, которая имеет математическое ожидание, в следующем виде:

(где h — среднее значение для распределения p (T), R — наибольшее значение r).

Используя уравнение z мы предполагаем, что Марковская цепь сошлась. В следствие ее сходимости мы получили инвариантное распределение, для которого мы можем посчитать математическое ожидание:

Если уравнение z выполняется, то мы можем констатировать эргодичность цепи P. E. Rossi, G. M. Allenby and R. McCulloch Bayesian Statistics and Marketing 2005 John Wiley & Sons, Ltd. ISBN: 0−470−86 367−6 P 50−51. Тем не мене, уравнение z не всегда поможет нам проверить цепь на эргодичность, поскольку иногда необходимо выполнять проверку, когда r довольно мало. Для проверки эргодичности цепи с малым числом итераций используется уравнение детального баланса.

Если уравнение z выполняется, мы также можем диагностировать эргодичность цепи P. E. Rossi, G. M. Allenby and R. McCulloch Bayesian Statistics and Marketing 2005 John Wiley & Sons, Ltd. ISBN: 0−470−86 367−6 p.88.

Процедура проверки сходимости и период «burn-in».

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

Рисунок 1. Средние значения частичных полезностей в по всем респондентам.

Для данного иерархического байесовского анализа была сгенерирована Марковская цепь в 20 тысяч итераций. Значения параметров сохранялись каждую 20-ю итерацию, поэтому на диаграмме всего тысяча средних значений Бета. График демонстрирует, что цепь достигла сходимости на 6-ой тысяче итераций. Остальные итерации с 1 по 6000 были сгенерированы для сходимости. Их можно удалить или «сжечь». В литературе период до сходимости называется периодом «burn in» P. E. Rossi, G. M. Allenby and R. McCulloch Bayesian Statistics and Marketing 2005 John Wiley & Sons, Ltd. ISBN: 0−470−86 367−6 P 51.

Алгоритмы для получения новых значений Схема Метрополиса-Гастингса Предположим, что нам необходимо сгенерировать выборку из распределения p0(T). Рассмотрим шаг генерации значений по схеме Метрополиса-Гастингса. Пусть на шаге-n сгенерирована конфигурация Tr, тогда на шаге r+1 сначала генерируется конфигурация T* из исходного распределения r (T|Tr). После того, как значение сгенерировано, вычисляется величина g:

Если g больше или равно 1, то значение T* принимается в качестве точки Tr для следующего шага с вероятностью g. В противном случае,.

Если, то цепь сохраняет инвариантность, поскольку значение Tr не изменяется. Для второго случая, когда необходимо провести проверку на эргодичность. В данном случае r может быть мало, поэтому для проверки мы будем использовать уравнение детального баланса z. Если цепь прошла проверку на эргодичность, то мы имеем право принять данную точку c вероятностью пропорциональной увеличению плотности. Если цепь не прошла проверку, то мы принимаем данную точку с вероятностью меньшей, чем данная плотность.

Схема Гиббса.

Сэмплер Гиббса позволяет сгенерировать выборку из многомерного распределения p (T), где T={}. Задача генерации выборки из многомерного распределения является задачей высшей степени сложности. Тем не менее, метод Гиббса позволяет обойти сложности и превращает процедуру генерации в более оптимальный процесс. Рассмотрим итерацию по схеме Гиббса. Допустим, на итерации r сгенерирован набор новых значений параметров Tr={ }. Тогда на следующей итерации алгоритм работает следующим образом:

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

Подводя итог по теоретическим нюансам Марковских цепей Монте Карло, в основе алгоритма MCMC лежит Марковская цепь. Идея цепи Маркова заключается в последовательной генерации новых распределений параметров при условии распределений этих же параметров на предыдущей итерации. Для генерации параметров из многомерного распределения используется сэмплер Гиббса. В процессе генерации параметров мы можем столкнуться с тем, что распределения новых параметров сомнительны. Для того, чтобы принять или отвергнуть распределения новых параметров следует воспользоваться алгоритмом Метрополис Гастингса. В течение всей процедуры MCMC необходимо следить за тем, чтобы цепь сохраняла эргодичность. Невыполнение условия эргодичности цепи приводит к генерации неинвариантных распределений параметров, на основании которых мы не можем посчитать корректные математические ожидания для выборки. Впоследствии мы не сможем получить интерпретируемые индивидуальные полезности для conjoint’a. Напротив, если наша цепь сошлась и соблюдает критерий эргодичности, то мы получаем корректные индивидуальные полезности. Чтобы окончательно убедиться в том, что мы получили корректные значения параметров, необходимо свериться с графиком средних частичных полезностей и отсечь период «burn in».