Интеграция Монте-Карло

Иллюстрация интеграции Монте-Карло. В этом примере областью D является внутренний круг, а областью E — квадрат. Поскольку площадь квадрата (4) можно легко вычислить, площадь круга (π*1,0 ²) можно оценить по отношению (0,8) точек внутри круга (40) к общему количеству точек (50), что дает приближенное значение площади круга 4 * 0,8 = 3,2 ≈ π.

В математике интегрирование Монте-Карло представляет собой метод численного интегрирования с использованием случайных чисел . Это особый метод Монте-Карло , который численно вычисляет определенный интеграл . В то время как другие алгоритмы обычно оценивают подынтегральную функцию в регулярной сетке, ^[1] Монте-Карло случайным образом выбирает точки, в которых вычисляется подынтегральная функция. ^[2] Этот метод особенно полезен для многомерных интегралов. ^[3]

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

Обзор

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

Проблема интегрирования Монте-Карло заключается в вычислении многомерного определенного интеграла. $I=\int _{\Omega }f({\overline {\mathbf {x} }})\,d{\overline {\mathbf {x} }}$ где Ω, подмножество $\mathbb {R} ^{m}$ , имеет объем $V=\int _{\Omega }d{\overline {\mathbf {x} }}$

Наивный подход Монте-Карло заключается в равномерной выборке точек на Ω: ^[4] учитывая N однородных выборок, ${\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in \Omega ,$

Меня можно приблизительно оценить $I\approx Q_{N}\equiv V{\frac {1}{N}}\sum _{i=1}^{N}f({\overline {\mathbf {x} }}_{i})=V\langle f\rangle .$

Это происходит потому, что закон больших чисел гарантирует, что $\lim _{N\to \infty }Q_{N}=I.$

Учитывая оценку I по Q _N , планки ошибок Q _N можно оценить по выборочной дисперсии, используя несмещенную оценку дисперсии .

$\mathrm {Var} (f)\equiv \sigma _{N}^{2}={\frac {1}{N-1}}\sum _{i=1}^{N}\left(f({\overline {\mathbf {x} }}_{i})-\langle f\rangle \right)^{2}.$ что приводит к $\mathrm {Var} (Q_{N})={\frac {V^{2}}{N^{2}}}\sum _{i=1}^{N}\mathrm {Var} (f)=V^{2}{\frac {\mathrm {Var} (f)}{N}}=V^{2}{\frac {\sigma _{N}^{2}}{N}}.$

Пока последовательность $\left\{\sigma _{1}^{2},\sigma _{2}^{2},\sigma _{3}^{2},\ldots \right\}$ ограничена, эта дисперсия асимптотически убывает до нуля как 1/ N . Таким образом, оценка ошибки Q _N равна $\delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sigma _{N}}{\sqrt {N}}},$ который уменьшается как ${\tfrac {1}{\sqrt {N}}}$ . Это стандартная ошибка среднего значения, умноженного на $V$ . Этот результат не зависит от количества измерений интеграла, что является обещанным преимуществом интегрирования Монте-Карло по сравнению с большинством детерминированных методов, которые экспоненциально зависят от размерности. ^[5] Важно отметить, что, в отличие от детерминистических методов, оценка ошибки не является строгой границей ошибки; случайная выборка может не раскрыть все важные особенности подынтегральной функции, что может привести к недооценке ошибки.

Хотя наивный метод Монте-Карло работает для простых примеров, улучшение по сравнению с детерминированными алгоритмами может быть достигнуто только с помощью алгоритмов, которые используют выборочные распределения для конкретной задачи. При соответствующем распределении выборки можно использовать тот факт, что почти все подынтегральные выражения более высокой размерности очень локализованы и только небольшое подпространство вносит заметный вклад в интеграл. ^[6] Большая часть литературы по методу Монте-Карло посвящена разработке стратегий улучшения оценок ошибок. В частности, двумя примерами таких методов являются стратифицированная выборка (разделение региона на поддомены) и выборка по важности (выборка из неравномерных распределений).

Пример

Парадигматическим примером интеграции Монте-Карло является оценка π. Рассмотрим функцию $H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}$ и множество Ω = [−1,1] × [−1,1] с V = 4. Заметим, что $I_{\pi }=\int _{\Omega }H(x,y)dxdy=\pi .$

Таким образом, грубый способ вычисления значения π с помощью интегрирования Монте-Карло состоит в том, чтобы выбрать N случайных чисел на Ω и вычислить $Q_{N}=4{\frac {1}{N}}\sum _{i=1}^{N}H(x_{i},y_{i})$

На рисунке справа относительная погрешность ${\tfrac {Q_{N}-\pi }{\pi }}$ измеряется как функция N , что подтверждает ${\tfrac {1}{\sqrt {N}}}$ .

Пример С

Имейте в виду, что следует использовать настоящий генератор случайных чисел.

int i, throws = 99999, insideCircle = 0;
double randX, randY, pi;

srand(time(NULL));

for (i = 0; i < throws; ++i) {
  randX = rand() / (double) RAND_MAX;
  randY = rand() / (double) RAND_MAX;
  if (randX * randX + randY * randY < 1) ++insideCircle;
}

pi = 4.0 * insideCircle / throws;

Пример Python

Сделано на Python .

from numpy import random

throws = 2000
inside_circle = 0
i = 0
radius = 1
while i < throws:
    # Choose random X and Y centered around 0,0
    x = random.uniform(-radius, radius)
    y = random.uniform(-radius, radius)
    # If the point is inside circle, increase variable
    if x**2 + y**2 <= radius**2:
        inside_circle += 1
    i += 1

# Calculate area and print; should be closer to Pi with increasing number of throws
area = (((2 * radius) ** 2) * inside_circle) / throws
print(area)

Пример Вольфрама Математики

Код ниже описывает процесс интеграции функции $f(x)={\frac {1}{1+\sinh(2x)\log(x)^{2}}}$ от $0.8<x<3$ используя метод Монте-Карло в системе Mathematica :

func[x_] := 1/(1 + Sinh[2*x]*(Log[x])^2);

(*Sample from truncated normal distribution to speed up convergence*)
Distrib[x_, average_, var_] :=   PDF[NormalDistribution[average, var], 1.1*x - 0.1];
n = 10;
RV = RandomVariate[TruncatedDistribution[{0.8, 3}, NormalDistribution[1, 0.399]], n];
Int = 1/n Total[func[RV]/Distrib[RV, 1, 0.399]]*Integrate[Distrib[x, 1, 0.399], {x, 0.8, 3}]

NIntegrate[func[x], {x, 0.8, 3}] (*Compare with real answer*)

Рекурсивная стратифицированная выборка

Иллюстрация рекурсивной стратифицированной выборки. В этом примере функция: $f(x,y)={\begin{cases}1&x^{2}+y^{2}<1\\0&x^{2}+y^{2}\geq 1\end{cases}}$
из приведенной выше иллюстрации, было проинтегрировано в пределах единичного квадрата с использованием предложенного алгоритма. Точки выборки были записаны и нанесены на график. Четко стратифицированный алгоритм выборки концентрирует точки в областях, где изменение функции является наибольшим.

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

Обычная стратегия «разделения на два» не работает для многомерных измерений, поскольку количество подтомов растет слишком быстро, чтобы их можно было отслеживать. Вместо этого оценивают, по какому измерению подразделение должно принести наибольшие дивиденды, и делят объем только по этому измерению.

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

Популярная программа MISER реализует аналогичный алгоритм.

МИЗЕР Монте-Карло

Алгоритм MISER основан на рекурсивной стратифицированной выборке . Этот метод направлен на уменьшение общей ошибки интегрирования за счет концентрации точек интегрирования в областях наибольшей дисперсии. ^[7]

Идея стратифицированной выборки начинается с наблюдения, что для двух непересекающихся областей a и b с оценками интеграла по методу Монте-Карло $E_{a}(f)$ и $E_{b}(f)$ и отклонения $\sigma _{a}^{2}(f)$ и $\sigma _{b}^{2}(f)$ , дисперсия Var( f ) комбинированной оценки $E(f)={\tfrac {1}{2}}\left(E_{a}(f)+E_{b}(f)\right)$ дается, $\mathrm {Var} (f)={\frac {\sigma _{a}^{2}(f)}{4N_{a}}}+{\frac {\sigma _{b}^{2}(f)}{4N_{b}}}$

Можно показать, что эта дисперсия минимизируется путем распределения точек таким образом, что: ${\frac {N_{a}}{N_{a}+N_{b}}}={\frac {\sigma _{a}}{\sigma _{a}+\sigma _{b}}}$

Следовательно, наименьшая оценка ошибки получается путем распределения точек выборки пропорционально стандартному отклонению функции в каждом субрегионе.

Алгоритм MISER делит область интегрирования пополам по одной оси координат, чтобы на каждом шаге получить две подобласти. Направление выбирается путем изучения всех d возможных делений пополам и выбора того, которое минимизирует совокупную дисперсию двух субрегионов. Дисперсия в субрегионах оценивается путем выборки с долей общего количества точек, доступных на текущем шаге. Затем та же процедура рекурсивно повторяется для каждого из двух полупространств наилучшего биссектрисы. Остальные точки выборки распределяются по субрегионам с использованием формулы для N _a и N _b . Это рекурсивное распределение точек интегрирования продолжается до заданной пользователем глубины, где каждый субрегион интегрируется с использованием простой оценки Монте-Карло. Эти отдельные значения и их оценки ошибок затем объединяются, чтобы получить общий результат и оценку его ошибки.

Выборка по важности

Существует множество алгоритмов выборки по важности, например:

Алгоритм выборки по важности

Выборка по важности представляет собой очень важный инструмент для интеграции Монте-Карло. ^[3]^[8] Основным результатом выборки по важности для этого метода является то, что единая выборка ${\overline {\mathbf {x} }}$ является частным случаем более общего выбора, при котором образцы берутся из любого распределения $p({\overline {\mathbf {x} }})$ . Идея в том, что $p({\overline {\mathbf {x} }})$ может быть выбран для уменьшения дисперсии измерения Q _N .

Рассмотрим следующий пример, в котором хотелось бы численно проинтегрировать функцию Гаусса с центром в 0 и σ = 1 от −1000 до 1000. Естественно, если выборки формируются равномерно на интервале [−1000, 1000], только очень малая часть из них будет иметь значение для интеграла. Это можно улучшить, выбрав распределение, отличное от того, из которого выбираются выборки, например, путем выборки в соответствии с гауссовым распределением с центром в 0 и σ = 1. Конечно, «правильный» выбор сильно зависит от подынтегральной функции.

Формально, учитывая набор образцов, выбранных из распределения $p({\overline {\mathbf {x} }}):\qquad {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in V,$ оценка для I дается выражением ^[3] $Q_{N}\equiv {\frac {1}{N}}\sum _{i=1}^{N}{\frac {f({\overline {\mathbf {x} }}_{i})}{p({\overline {\mathbf {x} }}_{i})}}$

Интуитивно это означает, что если мы выберем конкретный образец в два раза больше, чем другие образцы, мы весим его вдвое меньше, чем другие образцы. Эта оценка, естественно, справедлива для равномерной выборки, в случае, когда $p({\overline {\mathbf {x} }})$ является постоянным.

Алгоритм Метрополиса – Гастингса – один из наиболее часто используемых алгоритмов для генерации ${\overline {\mathbf {x} }}$ от $p({\overline {\mathbf {x} }})$ , ^[3] тем самым обеспечивая эффективный способ вычисления интегралов.

Вегас Монте-Карло

Алгоритм VEGAS аппроксимирует точное распределение, выполняя несколько проходов по области интегрирования, что создает гистограмму функции f . Каждая гистограмма используется для определения распределения выборки для следующего прохода. Асимптотически эта процедура сходится к искомому распределению. ^[9] Чтобы избежать роста количества интервалов гистограммы, как K ^д, распределение вероятностей аппроксимируется сепарабельной функцией: $g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots$ так что необходимое количество бункеров составляет всего Kd . Это эквивалентно расположению пиков функции по проекциям подынтегральной функции на оси координат. Эффективность VEGAS зависит от справедливости этого предположения. Это наиболее эффективно, когда пики подынтегральной функции хорошо локализованы. Если подынтегральную функцию можно переписать в приблизительно разделимой форме, это повысит эффективность интеграции с VEGAS. VEGAS включает ряд дополнительных функций и сочетает в себе как стратифицированную выборку, так и выборку по важности. ^[9]

См. также

Примечания

^ Пресс и др. 2007 , гл. 4
^ Пресс и др. 2007 , гл. 7
^ Jump up to: ^а ^б ^с ^д Ньюман и Баркема 1999 , гл. 2
^ Ньюман и Баркема 1999 , гл. 1
^ Пресс и др. 2007 год
^ Маккей 2003 , стр. 284–292.
^ Press & Farrar 1990 , стр. 190–195.
^ Крозе, Таймре и Ботев, 2011 г.
^ Jump up to: ^а ^б Лепаж 1978 г.

Ссылки

Кафлиш, RE (1998). «Методы Монте-Карло и квази-Монте-Карло». Акта Нумерика . 7 :1–49. Бибкод : 1998AcNum...7....1C . дои : 10.1017/S0962492900002804 . S2CID 5708790 .
Вайнцирль, С. (2000). «Введение в методы Монте-Карло». arXiv : hep-ph/0006269 .
Пресс, WH; Фаррар, Г. Р. (1990). «Рекурсивная стратифицированная выборка для многомерной интеграции Монте-Карло» . Компьютеры в физике . 4 (2): 190. Бибкод : 1990ComPh...4..190P . дои : 10.1063/1.4822899 .
Лепаж, врач общей практики (1978). «Новый алгоритм адаптивной многомерной интеграции». Журнал вычислительной физики . 27 (2): 192–203. Бибкод : 1978JCoPh..27..192L . дои : 10.1016/0021-9991(78)90004-9 .
Лепаж, врач общей практики (1980). «VEGAS: Адаптивная программа многомерной интеграции». Корнеллский препринт CLNS 80-447 .
Хаммерсли, Дж. М.; Хэндскомб, округ Колумбия (1964). Методы Монте-Карло . Метуэн. ISBN 978-0-416-52340-9 .
Крозе, ДП ; Таймре, Т.; Ботев, З.И. (2011). Справочник по методам Монте-Карло . Джон Уайли и сыновья.
Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007). Численные рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN 978-0-521-88068-8 .
Маккей, Дэвид (2003). «Глава 4.4 Типичность и глава 29.1» (PDF) . Теория информации, логический вывод и алгоритмы обучения . Издательство Кембриджского университета. ISBN 978-0-521-64298-9 . МР 2012999 .
Ньюман, МЭЖ; Баркема, GT (1999). Методы Монте-Карло в статистической физике . Кларендон Пресс.
Роберт, CP; Казелла, Дж. (2004). Статистические методы Монте-Карло (2-е изд.). Спрингер. ISBN 978-1-4419-1939-7 .

Внешние ссылки

Café math: Интеграция Монте-Карло : статья в блоге, описывающая интеграцию Монте-Карло (принцип, гипотеза, доверительный интервал).
Boost.Math: Наивная интеграция Монте-Карло: документация для наивных подпрограмм Монте-Карло на C++.
Апплет Монте-Карло, применяемый в задачах статистической физики

[1] Пресс и др. 2007 , гл. 4

[2] Пресс и др. 2007 , гл. 7

[newman1999ch2-3] Jump up to: ^а ^б ^с ^д Ньюман и Баркема 1999 , гл. 2

[newman1999ch1-4] Ньюман и Баркема 1999 , гл. 1

[5] Пресс и др. 2007 год

[6] Маккей 2003 , стр. 284–292.

[7] Press & Farrar 1990 , стр. 190–195.

[kr11-8] Крозе, Таймре и Ботев, 2011 г.

[Lepage,_1978-9] Jump up to: ^а ^б Лепаж 1978 г.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]