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

Численная техника

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

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

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

Обзор

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

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

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

Я могу быть приблизительно определен $${\displaystyle I\approx Q_{N}\equiv V{\frac {1}{N}}\sum _{i=1}^{N}f({\overline {\mathbf {x} }}_{i})=V\langle f\rangle .}$$

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

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

$${\displaystyle \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}.}$$что приводит к $${\displaystyle \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}}.}$$

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

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

Пример

Относительная ошибка как функция количества образцов, показывающая масштабирование ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$

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

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

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

Пример С

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

int i , throws = 99999 , insideCircle = 0 ; double randX , randY , pi ;          srand ( время ( NULL ));для ( i = 0 ; i < throws ; ++ i ) { randX = rand () / ( double ) RAND_MAX ; randY = rand () / ( double ) RAND_MAX ; если ( randX * randX + randY * randY < 1 ) ++ insideCircle ; }                               пи = 4,0 * внутренний круг / броски ;      

Пример на Python

Сделано на Python .

из  numpy  импорт  случайныйthrows  =  2000 inside_circle  =  0 i  =  0 radius  =  1 while  i  <  throws : # Выбрать  случайные X и Y с центром около 0,0  x =  random.uniform  ( -radius , radius ) y = random.uniform ( -radius , radius ) # Если точка находится внутри круга, увеличить переменную if x ** 2 + y ** 2 < = radius ** 2 : inside_circle + = 1 i + = 1                  # Вычислить площадь и вывести ее на печать; она должна быть ближе к числу Пи с увеличением количества бросков площадь  =  ((( 2  *  радиус )  **  2 )  *  внутренний_круг )  /  броски печать ( площадь )

Пример Wolfram Mathematica

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

func [ x_ ] := 1 / ( 1 + Sinh [ 2 * x ] * ( Log [ x ]) ^ 2 );    (*Выборка из усеченного нормального распределения для ускорения сходимости*) 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 }] (*Сравните с реальным ответом*)    

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

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

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

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

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

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

MISER Монте-Карло

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

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

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

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

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

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

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

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

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

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

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

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

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

ВЕГАС Монте-Карло

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

Смотрите также

• Метод квази-Монте-Карло
• Вспомогательное поле Монте-Карло
• Метод Монте-Карло в статистической физике
• Метод Монте-Карло
• Сокращение дисперсии

Примечания

1. Пресс и др. 2007, Гл. 4
2. Пресс и др. 2007, Гл. 7
3. Ньюман и Баркема 1999, Гл. 2
4. Ньюман и Баркема 1999, Гл. 1
5. Пресс и др. 2007
6. Маккей 2003, стр. 284–292.
7. Пресс и Фаррар 1990, стр. 190–195.
8. Крез, Таймре и Ботев 2011
9. Лепаж 1978

Ссылки

• Кафлиш, Р.Э. (1998). «Монте-Карло и квази-Монте-Карло». Акта Нумерика . 7 : 1–49. Бибкод : 1998AcNum...7....1C. дои : 10.1017/S0962492900002804. S2CID  5708790.
• Вайнцирль, С. (2000). «Введение в методы Монте-Карло». arXiv : hep-ph/0006269 .
• Press, WH; Farrar, GR (1990). "Рекурсивная стратифицированная выборка для многомерной интеграции Монте-Карло". Компьютеры в физике . 4 (2): 190. Bibcode :1990ComPh...4..190P. doi : 10.1063/1.4822899 .
• Лепаж, ГП (1978). «Новый алгоритм адаптивной многомерной интеграции». Журнал вычислительной физики . 27 (2): 192–203. Bibcode : 1978JCoPh..27..192L. doi : 10.1016/0021-9991(78)90004-9.
• Лепаж, ГП (1980). «VEGAS: адаптивная многомерная программа интеграции». Препринт Корнелла CLNS 80-447 .
• Hammersley, JM; Handscomb, DC (1964). Методы Монте-Карло . Methuen. ISBN 978-0-416-52340-9.
• Крезе, Д.П .; Таймре, Т.; Ботев, З.И. (2011). Справочник по методам Монте-Карло . John Wiley & Sons.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). Numerical Recipes: The Art of Scientific Computing (3-е изд.). Нью-Йорк: Cambridge University Press. ISBN 978-0-521-88068-8.
• MacKay, David (2003). "глава 4.4 Типичность и глава 29.1" (PDF) . Теория информации, вывод и алгоритмы обучения. Cambridge University Press. ISBN 978-0-521-64298-9. МР  2012999.
• Newman, MEJ; Barkema, GT (1999). Методы Монте-Карло в статистической физике . Clarendon Press.
• Роберт, CP; Казелла, Дж. (2004). Статистические методы Монте-Карло (2-е изд.). Спрингер. ISBN 978-1-4419-1939-7.

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

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