Глава 1. Метод Монте-Карло

Алгоритм Метрополиса в вопросах интегрирования методом Монте-Карло

Выполнил

студент 3 курса, группы ПМИ-31

Булдаков С. А.

_______________

Научный руководитель

канд. пед. наук

Соколова А. Н.

________________

Киров

Оглавление

Введение. 3

Глава 1. Метод Монте-Карло. 4

1.1. Понятие метода. 4

1.2. Алгоритм Метрополиса. 5

Глава 2. Реализация алгоритма Метрополиса. 6

2.1. Программная реализация алгоритма. 6

2.3. Анализ результатов эксперимента. 6

Заключение. 8

Библиографический список. 9

Приложение. 10

Введение

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

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

Для достижения поставленной цели должны быть решены следующие задачи:

1. Изучить теоретические сведения об алгоритме Метрополиса.

2. Написать программу, реализующую метод Монте-Карло для численного интегрирования.

3. Провести анализ результатов тестирования программы и сформулировать выводы.

Глава 1. Метод Монте-Карло

1.1. Понятие метода

Рассмотрим систему из N частиц описываемой функцией гамильтониана , где задает все степени свободы одной частицы (например ). Вектор задает одно состояние системы. Множество состояний системы составляет доступное ей фазовое пространство Ω. Тогда среднее значение величины A, являющейся функцией состояния системы, дается интегралом

(1.1)

где p — функция распределения, а в знаменателе находится статистическая сумма (1.2)

Если система состоит из небольшого числа частиц и размерность пространства Ω мала, то интеграл (1.1) можно вычислить, используя обычные формулы для приближенного численного вычисления интегралов с заданной точностью. Однако при большом N, когда кратность интеграла становится большой, такой подход малопродуктивен, т.к. затраты на вычисления зависят экспоненциально от N.

Другой способ, носящий имя метода Монте Карло, основан на стохастическом переборе точек в фазовом пространстве с предпочтительной выборкой тех областей из Ω, которые дают существенный вклад в интеграл (1.1). Таким образом, в соответствии с функцией распределения p генерируется цепь состояний в фазовом пространстве.

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

Обычно генерируется марковская цепь, т.е. такая последовательность, в которой последующее состояние зависит только от настоящего состояния и не зависит «от прошлого». Математически это означает, что условная вероятность P(R_n|Rn−1, ..., R0) появления состояния Rn после последовательности Rn−1, ..., R0 равна вероятности P(Rn|Rn−1).

Примером Марковской цепи может служить [5].

Одномерное дискретное случайное блуждание — это случайный процесс с дискретным временем, имеющий вид

,

где

— начальное состояние;

;

случайные величины совместно независимы.

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

Целью алгоритма Метрополиса является создание такой последовательности состояний a1, a2, ..., am, что в пределе m → ∞ распределение становится больцмановским.

Более подробную информацию о цепях Маркова можно получить, например, в [15], [6]. Сформулируем теперь сам алгоритм Метрополиса. Вероятность P(b ← a) принятия пробного движения преводящего систему из состояния a в состояние b задается следующим образом:

,

где p(a), p(b) – вероятности состояний a и b соответственно.

Можно показать, условие детального баланса выполнено и распределение состояний a1, a2, ... в фазовом простанстве дается функцией p(R). Таким образом, переход осуществляется в то из состояний a или b, которое имеет большую вероятность.

Очередная итерация алгоритма начинается с состояния и заключается в следующем:

• выбрать x' по распределению ;
• вычислить

;

• С вероятностью a (1, если a>1) , иначе .

Суть в том, что мы переходим в новый центр распределения, если примем очередной шаг, и получается случайное блуждание, зависящее от распределения p: мы всегда принимаем шаг, если в эту сторону функция плотности увеличивается, и иногда отвергаем, если уменьшается (отвергаем не всегда, потому что нам надо уметь выходить из локальных максимумов и блуждать под всем графикомp). Замечу, кстати, что когда шаг отвергается, мы не просто отбрасываем новую точку, а повторяем второй раз то же самое x⁽t⁾ — таким образом высока вероятность повторять сэмплы вокруг локальных максимумов распределения p, что соответствует большей плотности (рис. 1).

Рис. 1. Блуждание на плоскости под графиком смеси двух двумерных гауссианов