CC BY-NC-ND
223
О методах Монте-Карло с цепями Маркова
Шведов A.C.
В работе рассматриваются два метода Монте-Карло с цепями Маркова, широко применяемые в эконометрических исследованиях. Это алгоритм Метрополиса и гиббсовский выбор. Приводится описание обоих методов. Методы Монте-Карло с цепями Маркова предназначены для симулирования наборов векторов, отвечающих многомерным распределениям вероятностей. В частности, эти методы применяются в байесовской статистике для исследования апостериорных распределений. Существенное значение имеет соблюдение условия инвариантности, доказательства, что это условие выполняется, приводятся для обоих методов. Для обоснования и изучения методов используется теория цепей Маркова с конечным числом состояний. На нескольких примерах исследуется точность рассматриваемых методов Монте-Карло с цепями Маркова. Эти примеры включают двумерное нормальное распределение с высокой корреляцией, двумерное экспоненциальное распределение, смесь двумерных нормальных распределений.
Ключевые слова: симулирование многомерных распределений; алгоритм Метрополиса; гиббсовский выбор.
1. Введение
В последние два десятилетия в эконометрическую науку прочно вошли методы Монте-Карло с цепями Маркова, так называемые MCMC (Markov chain Monte Carlo) методы. Большое значение для развития этого направления эконометрики имеют современные высокоскоростные компьютеры. Математическая модель любой системы должна быть оправданно сложной, чтобы включать в себя основные черты этой системы, но не должна быть излишне сложной. Как правильно определить границу между оправданной и излишней сложностью? И в естественных науках, и в науках о человеке и обществе эта граница при решении многих трудных задач смещается в сторону все более сложных математических моделей. Во многом это определяется увеличением вычислительных мощностей, используемых при анализе.
Но чем более сложная эконометрическая модель используется, тем труднее ее исследовать аналитическими методами и тем чаще приходится прибегать к чис-
Шведов A.C. - д. физ.-мат. н., профессор кафедры математической экономики и эконометрики Государственного университета - Высшей школы экономики, e-mail: [email protected]
Статья поступила в Редакцию в январе 2010 г.
ленным методам. Примеры интересных и важных эконометрических моделей, использование которых становится возможным благодаря привлечению методов Монте-Карло с цепями Маркова, можно найти, скажем, в книгах [5] и [9]. Так, в [5] приводится пример, относящийся к аграрной экономике, в [9] - пример из портфельной теории. Число подобных примеров достаточно велико.
При применении метода Монте-Карло требуется построить набор случайных чисел, отвечающих тому или иному распределению вероятностей. Для одномерных функций распределения существующие подходы к построению наборов случайных чисел излагаются, например, в [1, гл. 4]. Такую же задачу построения набора теперь уже не случайных чисел, а векторов необходимо решать и для многомерных функций распределения. Для решения этой задачи и используются методы Монте-Карло с цепями Маркова.
Многие применения методов Монте-Карло с цепями Маркова в эконометрике связаны с байесовским подходом, когда эти методы используются для изучения апостериорного распределения параметров модели. Именно с таких позиций чаще всего ведется изложение и в обзорах, посвященных методам Монте-Карло с цепями Маркова. Но суть этих методов никак не связана с байесовским подходом. Одна из целей настоящей работы состоит в том, чтобы дать обзор соответствующих результатов в чистом виде, не увязывая с их возможным применением в байесовской эконометрике. Непростые в идейном отношении методы заслуживают такого аккуратного рассмотрения.
Кроме того, работа содержит некоторые новые соображения, касающиеся обоснования методов Монте-Карло с цепями Маркова.
Еще одна цель данной работы связана с достаточно длительной дискуссией о так называемой диагностике сходимости при применении методов Монте-Карло с цепями Маркова (см., например, [9]). Дискуссия вызвана практически важными вопросами о том, какой длины брать цепь, или, что лучше - использовать одну длинную цепь или несколько более коротких. В настоящей работе предлагается ответ с позиций, принятых в вычислительной математике. Бессмысленно говорить о какой-то диагностике сходимости, например, при приближенном вычислении интеграла с помощью квадратурной формулы. Чем больше используется узлов, тем точнее интеграл можно вычислить. Несколько подробно разобранных в данной статье примеров показывают, что в методах Монте-Карло с цепями Маркова ситуация совершенно такая же. Чем длиннее цепь, тем выше точность. На этих же примерах, для которых известно либо аналитическое решение, либо решение, посчитанное с большой точностью другим численным методом, исследуется, какой длины необходима цепь, чтобы обеспечить ту или иную точность ответа.
Пусть /(х) - совместная функция плотности компонент т -мерного случайного вектора. Рассматривается задача построения набора т -мерных векторов, отвечающих распределению вероятностей с функцией плотности /. Обозначим через В
множество точек х е Rm таких, что /(х) > 0 при х е В и /(х) = 0 при х П В.
Рассмотрим последовательность т -мерных случайных векторов Х0,Х1,...,Хп,..., которые составляют цепь Маркова. Начальные сведения о цепях Маркова можно найти, например, в [1, гл. 6]. Для этой цепи Маркова и для борелевского множества
А н Rm обозначим через Р("\х, А) вероятность перехода из состояния х в множество А ровно за п шагов. (Условие, что множество А является борелевским, накла-
дывается для того, чтобы по этому множеству можно было интегрировать. Параллелепипеды, шары, цилиндрические, конические множества - далеко не полное перечисление видов борелевских множеств. Примеры неборелевских множеств весьма искусственны и имеют отдаленное отношение к экономическим исследованиям. Но
все же брать в качестве А произвольное множество из Rm нельзя.)
Если для любого борелевского множества А н Rm при любом или почти любом х
то, реализовав цепь Маркова, можно получить набор т -мерных векторов, отвечающих распределению вероятностей с функцией плотности f. Обычное для одномерных распределений требование независимости случайных чисел, в данном случае, не накладывается. Об условиях, при которых имеет место (1), см., например, статью [2] и, особенно, теорему 6 из указанной статьи. Разумеется, для выполнения условия (1) цепь Маркова Х0,Х1,...,Хп,... должна каким-то образом строиться по функции f.
Методы Монте-Карло с цепями Маркова применяются при работе не только с непрерывными, но и с дискретными и смешанными многомерными распределениями. В качестве примера здесь можно назвать статью [11]. Но в данной работе мы ограничиваемся непрерывными многомерными распределениями. В разделе 2 настоящей работы рассматривается алгоритм Метрополиса, в разделе 3 - гиббсовский выбор.
Этот алгоритм построения набора т -мерных векторов, отвечающих распределению вероятностей с функцией плотности f, был предложен (в несколько меньшей общности, чем излагается ниже) в работе [10]. Затем он изучался и развивался в работах [3], [8], [12], где и был доведен до современного уровня общности. Название алгоритма по фамилии только первого из авторов работы [10] является установившимся, хотя, видимо, не совсем справедливым.
Цепь Маркова Х0,Х1,...,Хп,... ищется среди тех, для которых вероятности перехода за один шаг для любого х е В и для любого борелевского множества А н Rm имеют вид
Здесь 1А (х) - индикатор множества А, т.е. 1А (х) = 1, если х е А, 1А (х) = 0, если х П А. Очевидно, что
(1)
2. Алгоритм Метрополиса
поскольку Р(1)( х, Rm) = 1.
Функция к(х, у) строится так, чтобы выполнялось условие реверсивности:
(4) /(х) к(х, у) = /(у) к(у, х)
для любых х,у е R'n.
Покажем, что из выполнения условия реверсивности (4) следует выполнение условия инвариантности:
(5) Т /(у) dy = Г тРх,А) /(х) dx
Л А
для любого борелевского множества А н R'n.
Действительно, используя (4) при втором переходе и (3) - при четвертом, получаем
Т тРт(х, А) /(х) dx = Т т/(х) g(х) 1А (х) dx + Т т Т /(х) к(х, у) dydx =
= Т / (х) g (х) йх + Т Т / (у) к(у, х) dy йх =
¿А ¿RmiЛ
= Т / (х) g (х) йх + Т Т / (у) к(у, х) йх йу = 2а
= ТА /(х) g(х) йх + £ /(у) (1 - g(у)) йу = £ /(у) йу, т.е. имеет место (5).
Как показано в [2], именно условие инвариантности (5) является основным для выполнения (1).
«Рычагом управления» в алгоритме Метрополиса является неотрицательная функция с(х, у), которую следует выбрать так, чтобы
(6) Т тс(х, у) йу
2Rm
=1
для любого х е Rm. Далее принимается
а( х, у) =
тщИ, /(у с(у х I при /(х) с(х,у) > 0 1 '/(х) с(х,у) * -м ^
1 при /(х) с(х, у) = 0.
Существенно отметить, что если х е В, у П В и с(х,у) > 0, то а(х,у) = 0. Проверим, что функция
(7) к(х,у) = с(х,у) а(х,у)
удовлетворяет условию реверсивности (4). Выберем точки х и у.
Пусть f (х) с(х, у) = 0. В этом случае левая часть в (4) равна нулю. Если f (У) с(у, х) = 0, то и правая часть в (4) равна нулю. Если f (у) с(у, х) > 0, то а (у, х) = 0, и вновь правая часть в (4) равна нулю. Таким образом, при f (х) с(х, у) = 0 условие 4 выполняется.
Пусть f (х) с(х,у) > 0. Если f (х) с(х,у) < f (у)с(у, х), то а(х,у) = 1,
а( у, х)=т^хА.
У (у) с( у, х)
Поэтому f(х) Цх,у) = f(х) с(х,у) = /(у)с(у,х)а(у,х) = f(y)Цу,х), и условие (4) выполняется. Если f (х) с(х, у) > f (у) с(у, х), то а (у, х) = 1,
а(х,у) = У(у)с(у,х).
/ (х) с( х, у)
Тогда У(х) ^х, у) = /(х) с(х, у) а(х, у) = У(у) с(у, х) = У(у) с(у, х)а(у, х) = У(у) ^у, х).
Показано, что функция ^х, у), задаваемая формулой (7), удовлетворяет условию реверсивности (4). Функция g (х) строится по формуле (3), заметим, что эта функция может принимать значения лишь между нулем и единицей.
Реализацию цепи Маркова Х0,Х1,...,Хп,... обозначим х(0), х(1),..., х(п),... Пусть как-то выбран вектор х(0) е В. В алгоритме Метрополиса вектор х("+Г) строится по
вектору х(-п) следующим образом.
Шаг 1. Симулируется т -мерный вектор у из распределения с плотностью
с(х(п), у).
Шаг 2. Симулируется число и из равномерного на отрезке [0,1] распределения. Шаг 3. При и < а(х("\у) принимается х(~п+1'> = у; при и > а(х("\у) принимается
х(п+1) = х(п).
Не при любом выборе функции с(х,у), удовлетворяющей условию (6), построенный набор т -мерных векторов х(0),х('Г),...,х(п\... отвечает распределению вероятностей с функцией плотности У. Например, если с(х, у) = 0 при любых х, у е В, то очевидно, что это не так. Приводимые ниже рассуждения помогают понять, каким условиям должна удовлетворять функция с(х, у).
В статье [2] доказательство того, что имеет место (1), - это сложная теорема. Нет сомнений, что должно существовать более простое, пусть и более грубое обоснование алгоритма Метрополиса. Какие соображения привели авторов работ [3, 10] к их алгоритму? Далее приводится интуитивное объяснение алгоритма Метрополиса. Чтобы понять это интуитивное объяснение, достаточно знать теорию цепей Маркова с конечным числом состояний. А для чтения статьи [2] необходимо знание существенно более сложной общей теории цепей Маркова.
Будем считать, что множество В ограничено. Предположим также, что это множество можно разбить на т -мерные кубы А1,...,Ам, длина стороны каждого из кубов е - малое положительное число. Положим
(8) лг =Т /(х)йх, г = 1,...,М.
А
Пусть у^,...,у(м) - центры кубов А1,...,Ам. Будем считать, что х(0) - это одна из точек у^,...,у{(М\ Далее предположим, что любая симулированная на шаге 1
точка у, если она оказалась принадлежащей кубу Аг, заменяется на у(~г\ Это приводит к новой цепи Маркова, имеющей лишь конечное число состояний. Эти состояния можно занумеровать числами 1,..., М.
Покажем, что для этой новой цепи Маркова с конечным числом состояний вероятности (8) приближенно являются инвариантными. (Определение инвариантных вероятностей для цепи Маркова с конечным числом состояний дается, например, в [1, гл. 6].)
При х('п) = у^г) в цепи Маркова с конечным числом состояний новое значение у1-*) получается на шаге 1 с вероятностью, приближенно равной етс(у (~г \ у(х)). Это значение оставляется с вероятностью а(у(г \ ). Поэтому с учетом (7) вероятность получить на шаге 3 новое значение у() приближенно равна
етс(у),у«) а(у(г),у«)= етк(у(г),у«).
Если новое значение не принимается, то цепь остается в состоянии у(~г\ Это происходит с вероятностью, приближенно равной
1 -е>тк(у(г ^ ум).
Заменяя сумму на интеграл, получаем другое приближенное выражение для этой вероятности:
1 -]>(у(г ^ у )йу.
В соответствии с (3) последнее выражение равно g (у (~г)).
Таким образом, переходные вероятности рга, т.е. вероятности Р(1 (у(~г\ {у(х)}), для цепи Маркова с конечным числом состояний приближенно могут быть определены по следующей формуле:
= íg (у(г))+ втк(у(г\у ^ ) при г = * Рг" \етк(у(г),у«) при г №
Чтобы найти уравнения, которым приближенно удовлетворяют вероятности
(8), вернемся к рассмотрению исходной цепи Маркова и воспользуемся формулой (5) при А = А5 (свойство инвариантности для исходной цепи Маркова установлено!). При А = А5 левая часть в формуле (5) равна я5. Для правой части формулы (5) с
учетом (2) и аппроксимации лг = ет /(у(~г)) может быть записана следующая цепочка приближенных равенств:
т (х) 1А„ (х) + ТА к(х, У) ёуЦ /(х) ёх = = Т g (х) / (х) ёх +Т Т И( х, у) / (х) ёхёу =
= ет g(/*)/(/*)+ Т ^^Уе /(у(Г))ёу =
* г
=g (у « к+ет *(у(г), у(х) Я=емп рр.
Итак, приближенно, для вероятностей лг получена следующая система линейных уравнений:
^Ы Ы
рг =1, = е ,Ршрг пРи ^ = ъ..^Ы.
г=1
Если все переходные вероятности рг5 для цепи Маркова с конечным числом состояний положительны, то системе уравнений (10) удовлетворяют инвариантные вероятности такие, что при любых г и 5
Р(п) (у(г), {уМ }) ® Р при п ® ю
(см., например, [1, гл. 6, теорема 6.1]). При условии положительности всех переходных вероятностей рг5 система линейных уравнений (10) имеет единственное решение. И если вероятности (8) этой системе линейных уравнений удовлетворяют, то они являются инвариантными для цепи Маркова с конечным числом состояний.
Это и требуется, но следует помнить, что соотношения, о которых говорится выше, имеют место не в точности, а приближенно.
Если обоснование алгоритма не является безупречно строгим (как, например, приведенное выше объяснение алгоритма Метрополиса), то на вопрос о пригодности или непригодности алгоритма, или насколько алгоритм хорош, ответ могут дать расчеты. В этой области математической науки велика роль вычислительного опыта.
Положительность всех переходных вероятностей рг5, задаваемых формулой
(9), обеспечивается условием с(х, у) > 0 при любых х, у о В. Однако практически функцию с(х,у) можно брать с некоторым «запасом», чтобы условие с(х,у) > 0 выполнялось и при каких-то у П В, если это помогает в проведении шага 1.
Можно указать на аналогию между парой «исходная цепь Маркова и аппроксимирующая ее цепь Маркова с конечным числом состояний» и парой «дифферен-
циальное уравнение и разностная схема». Однако эта аналогия не является полной. В данном случае для расчета используется исходная цепь Маркова, а не цепь Маркова с конечным числом состояний.
Пример 1. Пусть / - функция плотности двумерного нормального распределения
(11)
7 (х)=ехр(~ 2(х - м)'Е-1 (х - т)
где
т =
Е =
а1
_2
а
2 Ш
(см., например, [1, гл. 1, § 6]). Штрих означает транспонирование.
Разумеется, чтобы построить набор векторов, отвечающих двумерному нормальному распределению, нет необходимости использовать метод Монте-Карло с
цепью Маркова. Достаточно построить набор векторов г(0),г(1),...,г("\..., отвечающих двумерному нормальному распределению с нулевым средним и единичной ковариационной матрицей. (Задача построения таких векторов сводится к задаче построения стандартных нормальных случайных чисел. Различные подходы к построению нормальных случайных чисел рассматриваются, например, в [1, гл. 4, § 2].) Затем, взяв матрицу
L =
V1 - р2°1
Р°1
0
2 Ш
для которой выполняется соотношение Е = Ь' Ь, надо положить
(12)
х(п) = т + Ь'г(п), п = 1,2,...
Но мы используем функцию плотности (11), чтобы проиллюстрировать алгоритм Метрополиса.
При каждом а > 0 на плоскости Я2 с координатами х1, х2 можно рассмотреть эллипс рассеивания
(х - /) Е-1 (х - т) = а2,
где х = (х1, х2) . Тогда вероятность того, что двумерный нормальный вектор с функцией плотности (11) принимает значение внутри этого эллипса, равна
(13) 1 - ехр(- 0,5 а2)
(см., например, [1, гл. 1, § 6]).
В расчете принимаются значения ¡и1 = 0, т2 = 0, = 0,8, о2 = 1,2, р = 0,9. Вероятность р1 рассчитывается при каждом а по формуле (13).
Во-первых, при различных N по формуле (12) строится набор точек Через р2 обозначается доля точек из этого набора, попавших внутрь эллипса рассеивания.
Во-вторых, при помощи алгоритма Метрополиса реализуется цепь Маркова. Используется функция с(х, у) = /0(у - х), где /0 - функция плотности двумерного нормального распределения с нулевым средним и с ковариационной матрицей
0,36 0 Ц 0 0,364.
Принимается начальное значение х1-0-1 = 0, и первые 500 построенных точек игнорируются. Таким образом, рассматривается набор
точек х(500+1),..., х(500Через р3 обозначается доля точек из этого набора, попавших внутрь эллипса рассеивания. Результаты расчетов представлены в табл. 1.
Таблица 1.
Соответствие долей точек, попавших внутрь эллипсов рассеивания, построенных различными методами (р2 — с использованием аналитической формулы (12), р3 — с использованием алгоритма Метрополиса), и вероятностей р1 попадания точки внутрь эллипса рассеивания, рассчитанных по формуле (13),
при различных N
а Р1 Р2 Р1 - Р2 Р3 Р1 - Р3
а) N = 500
0,5 0,11750 0,09000 0,02750 0,12000 -0,00250
1,0 0,39347 0,37000 0,02347 0,45400 -0,06053
1,5 0,67535 0,66000 0,01535 0,73000 -0,05465
2,0 0,86466 0,87200 -0,00734 0,85200 0,01266
Ь) N = 5000
0,5 0,11750 0,10900 0,00850 0,13640 -0,01890
1,0 0,39347 0,38060 0,01287 0,42780 -0,03433
1,5 0,67535 0,66960 0,00575 0,72160 -0,04625
2,0 0,86466 0,86040 0,00426 0,87600 -0,01134
с) N = 50000
0,5 0,11750 0,11708 0,00042 0,12728 -0,00978
1,0 0,39347 0,39000 0,00347 0,39866 -0,00519
1,5 0,67535 0,67392 0,00143 0,68304 -0,00769
2,0 0,86466 0,86232 0,00234 0,86420 0,00046
Окончание табл. 1.
а Р1 Р2 Р1 - Р2 Р3 Р1 - Рз
d) N = 500000
0,5 0,11750 0,11722 0,00028 0,11779 -0,00029
1,0 0,39347 0,39268 0,00079 0,39470 -0,00123
1,5 0,67535 0,67561 -0,00027 0,67827 -0,00292
2,0 0,86466 0,86448 0,00018 0,86646 -0,00179
Из табл. 1 видно, что, как и следовало ожидать, построенная с помощью аналитической формулы (12) величина р2 лучше аппроксимирует вероятность Р1, чем построенная с помощью численного алгоритма величина р3. Однако и для алгоритма Метрополиса аппроксимация существенно улучшается при увеличении N. Пример 2. Пусть функция плотности имеет вид
/ (х) = 0,5 / т( х) + 0,5/(2)( х),
где /(1)(х) и /(Т)(х) - функции плотности двумерного нормального распределения (см. (11)) соответственно со средними /л('Г), М(2 и с ковариационными матрицами Е(1), 1,(2\ В расчете принимаются значения ¡¿р =-0,3, /ир = 0,5, о1(1) = 0,7, <7^ = 11, рт = 0,8, тР = 0,8, =-0,2, а1(2) = 0,9, а^ = 1,0, р(2) =-0,1. Для построения набора векторов, отвечающих распределению вероятностей с функцией плотности /, используется алгоритм Метрополиса. Берется та же функция с(х, у), что и в примере 1, принимается начальное значение х1-0-1 = 0, и первые 500 построенных точек игнорируются.
Полное исследование того, насколько хорошо по набору точек х'
(2),
(500+1) х(500 +Ы)
можно судить о свойствах данного распределения, мы не проводим. Сравним лишь
Р1 = Т Т /(х1,х2)йх1йх2
и долю точек из рассматриваемого набора, попавших в квадрат [0,1] х [0,1]. Эту долю обозначим через р2. Применение кубатурной формулы с достаточно мелкой сеткой показывает, что р1 приближенно равно 0,10614. Полученные в расчете значения р2 при различных N приводятся в табл. 2.
Хотя самая высокая точность в данном расчете получена при N = 5000, расчеты с другими наборами исходных случайных чисел не дают повторения этого эффекта. Однако во всех расчетах видна примерно одна и та же тенденция к улучшению точности при увеличении N.
Таблица 2.
Соответствие долей точек р2, попавших внутрь квадрата [0,1] х [0,1], построенных с использованием алгоритма Метрополиса, и вероятности р1 попадания точки внутрь данного квадрата, рассчитанной по кубатурной формуле с большим числом узлов
N
Р2
Р1 - Р2
500 5000 50000 500000
0,08200 0,10600 0,10528 0,10640
0,02414 0,00014 0,00086 -0,00026
3. Гиббсовский выбор
Здесь описывается еще один алгоритм построения набора т -мерных векторов, отвечающих распределению вероятностей с функцией плотности /. Название «гиббсовский выбор» для данного алгоритма предложено в статье [7], поскольку используется он в этой работе для распределения, носящего имя выдающегося физика конца XIX - начала XX вв. Дж.У. Гиббса. Данное название закрепилось после того, как этот алгоритм был представлен в работе [6]. О более ранних работах, где используются сходные алгоритмы, см. [5, с. 20] или [9, с. 12].
Пусть х = (х1,...,хк), где хг е Ят, т1 +... + тк = т, и пусть множество В имеет
вид
В = В1 х ...х Вк,
где Вг с Ят - борелевское множество положительной меры.
При I = I,...,к рассмотрим маргинальную функцию плотности
/г хк ) = Т хг-Ъ Яг, хг+1,..., хк ) ¿Яг
2 в,
(при г = 1 отсутствуют аргументы х1,...,хг-1; при г = к отсутствуют аргументы хг+1,...,хк ). Фактически является функцией аргументов х1,...,хг-1, хг+1,...,хк.
Тогда при х е В можно определить условную функцию плотности gi (х) соотношением
/(х) = gi(х)/г(х), г = 1,...,к. Используя распространенное обозначение, можно сказать, что gi(х1,...,хк) - это условная плотность
хг \ х1,...,хг-1, хг+1,...,хк
(см., например, [1, гл. 1, § 3]). Будем считать, что gi (х) = 0 при х П В.
Цепь Маркова Х0,Х1,...,Хп,... ищется среди тех, для которых вероятности перехода за один шаг для любого х е В и для любого борелевского множества А н Rm имеют вид
Р (1)( х, А) = Т К( х, у) dy
Л А
(ср. (2)). Положим
(14) К( х у) = Х & (у^" У', х'+l,..., хк)
и проверим, что выполняется условие инвариантности (5). Имеем
/(У1,...У'-их',х'+1—хк)= 8'(У1,...У'-их',х'+1^хк) /'(У1,...У'-их',х'+1^хк) и аналогично
/(У1,...У'-иУ',х'+1^хк)= 8' (У1,...У'-иУ',х'+1^хк) /' (У1,...У'-иУ',х'+1^хк).
Тогда при условии, что аргументы функций /, и / являются точками множества В
/(У1,...У'-их',х'+1—хк) _ /(У1,...У'-иУ',х'+1—хк)
8;
(У1,...У'-и х', х'+1—хк) 8' (У1,...У'-иУ', х'+1^хкУ
Отсюда следует, что (15)
8' (У1,...У'-иУ', х'+1^.^хк)/(У1,...У'-и х', х'+1^.^хк) = = 8'(У1,...У'-и х', х'+l,..., хк) / (У1,...У'-и У', х'+1^ хк).
Данный прием исключения маргинальных плотностей с целью установить связь между условными и совместными плотностями, по-видимому, впервые был использован в [4, с. 195, 196]. Из формулы (15) нетрудно увидеть, что двум различным совместным функциям плотности /(1) и /(2 (при одном и том же множестве В) не может соответствовать один набор условных функций плотности 81,..., 8к.
Имеем
Т Р (1\ х, А) / (х) dx = Т Т / (х) К х, У) dydx = Т J (У) ау,
3 В ¡В1 А * А
где
J(у) = Т /(х) К(х, У) ах.
В
Из выражения (14) для ^х, у) находим
J(У) =Т gkСу^-Ук)Т gk-l(Уl,...,Ук-ихк)..Т gl(Уl,хк)f(xl,...,хк^...¿хк.
J Вк Вк-1 •>В1
Множитель gi (у1,...уг-1, уг, хг+1,..., хк ) намеренно не выносится из под знака интеграла по Вг. На основании (15) подынтегральная функция самого внутреннего интеграла равна
gl(x1, х2 хк ) 7 (У^ х2 хк ).
Теперь второй множитель можно вынести из под знака интеграла по В1 и воспользоваться тем, что
Т gl(xl, хк )ах1 = 1 •> В1
Тогда
J (У) =Т gk Су^- Ук )Т gk-1 Су^-^ Ук-Ъ хк )...Т g 2 (У^ У 2, хк ) f(Уl, хк ) ¿х2.. .¿хк.
J Вк Вк-1 •> В2
На основании (15) подынтегральная функция самого внутреннего интеграла равна
g 2 (Уl, х2 , хк ) 7 (Уи У 2 , ^З^ хк ).
Поэтому данный интеграл равен / (у1, у2, х3,..., хк ). Проведя аналогичные преобразования и для других подынтегральных функций, получаем 3(у) = /(у). Свойство инвариантности (5) установлено.
Как и в предыдущем разделе, реализацию цепи Маркова Х0,Х1,...,Хп,... обозначим х(0),х(),...,х(п\... Пусть как-то выбран вектор х(0) е В. При гиббсовском выборе вектор х(п+1) строится по вектору х(п) следующим образом.
Шаг 1. Симулируется т1 -мерный вектор у1 из распределения с функцией
плотности g1(y1, 4п),-, х^ ).
Шаг г. Симулируется тг1 -мерный вектор у из распределения с функцией
плотности gг (у^..., Уг-и Уг, x(i+)l,..., ^ ).
Шаг к. Симулируется тк -мерный вектор ук из распределения с функцией плотности gk УкУк ).
Принимается х("+г> = (у1,..., ук ).
Интуитивное объяснение гиббсовского выбора аналогично объяснению алгоритма Метрополиса.
Пример 3. Рассмотрим ту же функцию плотности двумерного нормального распределения (11), что и в примере 1. Условные плотности имеют вид
где
gl(xl, Х2 ) = "
о1с42л
ехр
g 2 (Хи Х2 )='
2&2
V (х2- т,с)2
а
2, с
422л
ехр
2 а
2— (х2 - М2, с )2
СТ1,с = ^л/1 - Р2 , °2, с = - Р
т\,с = т + р—(х2- т), тс = т + р—(х1- т)
а
(см., например, [1, гл. 1, § 6]).
Вектора х(\...,х(п\... строятся по следующим формулам. Пусть z1[\ 2'2>,...,,... - стандартные нормальные случайные числа. Тогда
_ (п+2)
^ 1
(16)
+2) = т + р ^ (х2п) - т)+^
х2п+Х) = т + р^2 (х2п+г> - т)+^2, с^2П+1.
Я
7 (п>
Интересно отметить, что отсюда следуют уравнения авторегрессии первого порядка
х(1п+1> = -Р2> + р2х2(п> + Б(п+л>, х(^+1> = т2(1 -р2> + р2х{п + е2п+1>.
Выражения для е!^" > и е^ > легко выводятся из (16).
Как и в примере 1, принимается начальное значение х(~а> = 0, и первые 500 построенных точек игнорируются. Таким образом, рассматривается набор точек . Через р4 обозначается доля точек из этого набора, попавших
х
(500+1>
..., х
(500+N >
внутрь эллипса рассеивания. Через р1 обозначается та же вероятность, что и в примере 1.
Результаты расчетов представлены в табл. 3.
Точность, с которой величина р4 аппроксимирует р1, улучшается с увеличением N. В расчетах, результаты которых приведены в данном примере и в примере 1, гиббсовский выбор дает более высокую точность, чем алгоритм Метрополиса. Гиббсовский выбор дает точность, которая близка к точности результатов, полученных при использовании аналитической формулы (12). Это можно объяснить тем, что в
1
с
1
2
данном примере гиббсовский выбор сводится к применению аналитических формул (16).
Таблица 3.
Соответствие долей точек р4, попавших внутрь эллипсов рассеивания, построенных с использованием гиббсовского выбора, и вероятностей р1 попадания точки внутрь эллипса рассеивания, рассчитанных по формуле (13), при различных N
а Р1 Р4 Р1 - Р4
а) N = 500
0,5 0,11750 0,12200 -0,00450
1,0 0,39347 0,39600 -0,00253
1,5 0,67535 0,69800 -0,02265
2,0 0,86466 0,88200 -0,01734
Ь) N = 5000
0,5 0,11750 0,12040 -0,00290
1,0 0,39347 0,39540 -0,00193
1,5 0,67535 0,67600 -0,00065
2,0 0,86466 0,87000 -0,00534
с) N = 50000
0,5 0,11750 0,11768 -0,00018
1,0 0,39347 0,38888 0,00459
1,5 0,67535 0,66866 0,00669
2,0 0,86466 0,86064 0,00402
d) N = 500000
0,5 0,11750 0,11741 0,00010
1,0 0,39347 0,39379 -0,00032
1,5 0,67535 0,67486 0,00048
2,0 0,86466 0,86445 0,00021
Пример 4. Пусть т = 2, множество В состоит из тех точек х = (х1, х2), для которых х1 > 0, х2 > 0. Совместная функция плотности имеет вид
/(х, х2) — СС ехр( х1х2 х х2),
где С - нормирующая константа. В этом примере гиббсовский выбор используется для расчета маргинальной функции плотности.
При т = 2 удобнее использовать другие обозначения для маргинальных функций плотности, чем те, которые применялись выше. Пусть
л(х1 )= Т Д^ х2) ¿х2^ 72 (х2 )= Т х2) ¿х1.
J 0 »0
Из вида / (х1, х2 ) следует, что
С
(17) /1(х1 ) =--ехр(- х1),
х1 +1
и аналогичная формула имеет место для /2 (х2).
Тогда, как нетрудно увидеть, условные функции плотности
gl (х1, х2 )= (х2 + 1)ехР(- (х2 + 1)х1 ), g2 (х1, х2 )= (х1 + 1)ехР(- (х1 + 1)х2 )
являются функциями плотности экспоненциального распределения.
Для определения значения функции /1 при некотором фиксированном х1 по
набору точек х(500+1),...,х(500+N), построенных при помощи гиббсовского выбора (как и в предыдущих примерах первые 500 точек игнорируются), используется следующее выражение:
(х )=1 еN * (х, х(5о°+п))
(18) Ш ) = ^ (х1, хГ^).
Поскольку
/(х1 )= Т gl(x1, х2) 72 (х2 ) ^
0
а правая часть в (18) - это приближенное выражение для правой части последнего равенства. Преимущества такого подхода к нахождению маргинальной функции плотности перед прямым подходом, когда расчеты производились бы по х^500х[500+N), объясняются, например, в [6].
В табл. 4 при различных N показаны значения
тах
х:ее
7(х1)- Л(х1)
где множество Q состоит из точек 0,04 • j при j = 1,...,100; функция /1 определяется по формуле (17), а функция 7 - по формуле (18); принимается начальное значение х(0) = (0,5,0,5).
Расчеты показывают, что точность аппроксимации улучшается с увеличением N.
Таблица 4.
Точность расчета маргинальной функции плотности при использовании гиббсовского выбора
N = 500 N = 5000 N = 50000 N = 500000
0,01875 0,00554 0,00360 0,00024
* * *
СПИСОК ЛИТЕРАТУРЫ
1. Шведов А.С. Теория вероятностей и математическая статистика - 2 (промежуточный уровень). М.: Изд. дом ГУ ВШЭ, 2007.
2. Athreya K.B., Doss H., Sethuraman J. On the Convergence of the Markov Chain Simulation Method // Annals of Statistics. 1996. 24. P. 69-100.
3. Barker A.A. Monte Carlo Calculations of the Radial Distribution Functions for a Proton-electron Plasma // Australian Journal of Physics. 1965. 18. P. 119-133.
4. Besag J. Spatial Interaction and Statistical Analysis of Lattice Systems // Journal of the Royal Statistical Societe. Series B. 1974. 36. P. 192-225.
5. Chen M.-H., Shao Q.-M, Ibragim J.G. Monte Carlo Methods in Bayesian Computation. Springer, 2000.
6. Gelfand A.E., Smith A.F.M. Sampling-based Approaches to Calculating Marginal Densities // Journal of the American Statistical Association. 1990. 85. P. 398-409.
7. Geman S., Geman D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images / / IEEE Transactions on Pattern Analysis and Machine Intelligence. 1984. 6. P. 721-741.
8. Hastings W.K Monte Carlo Sampling Methods Using Markov Chains and Their Applications // Biometrika. 1970. 57. P. 97-109.
9. Markov Chain Monte Carlo in Practice / W.R. Gilks, W.R.S. Richardson, D. Spiegelhalter (eds.) L.: Chapman and Hall, 1996.
10. Metropolis N, Rosenbluth A.W., Rosenbluth M.N., Teller A.H., Teller E. Equation of State Calculations by Fast Computing Machines // Journal of Chemical Physics. 1953. 21. P. 1087-1092.
11. Neal P., Rao T.S. MCMC for Integer-valued ARMA Processes // Journal of Time Series Analysis. 2007. 28. P. 92-110.
12. Peskun P.N. Optimum Monte Carlo Sampling Using Markov Chains // Biometrika. 1973. 60. P. 607-612.
