Меню Главная Случайная статья Настройки Метод Монте-Карло: различия между версиями Материал из https://ru.wikipedia.org Методы Монте-Карло (или алгоритмы численного статистического моделирования) – это алгоритмы реализации с помощью компьютера стохастической (вероятностной) модели того или иного объекта, с целью оценивания изучаемых интегральных параметров (средних значений нужных характеристик) на основе закона больших чисел. Методы используются для решения задач в различных областях физики, химии, математики, экономики, оптимизации, теории управления и др. Например для изучения случайных процессов. Содержание 1 История 1.1 Численное статистическое моделирование 1.2 Разработка теории и приложений алгоритмов 2 Основная часть 2.1 Общая схема метода Монте-Карло. Понятие монте-карловской оценки 2.2 Средняя погрешность и трудоемкость метода Монте-Карло 2.3 Методический пример: вычисление интеграла методом Монте-Карло 2.4 Компьютерное моделирование случайных величин, векторов и функций 3 Литература История Численное статистическое моделирование С развитием вычислительной техники возрастает роль компьютерных (численных) методов решения прикладных задач. Важное место в этом процессе занимают алгоритмы численного статистического моделирования или методы Монте-Карло. Особое место этих методов связано с простотой и естественностью их распараллеливания с целью эффективного применения современных многопроцессорных компьютерных систем. Под численным статистическим моделированием обычно понимают реализацию с помощью компьютера вероятностной модели того или иного объекта с целью оценивания изучаемых интегральных параметров (средних значений нужных характеристик) на основе закона больших чисел. Разработка теории и приложений алгоритмов Исторически интенсивное развитие теории и приложений метода Монте-Карло было связано с разработкой численных моделей ядерных процессов (при создании соответствующих военных и технических устройств — бомб, реакторов и т. п.) в СССР и США в 50-е годы XX века. Разработка теории метода связана с фундаментальными работами Дж. Неймана, С. Улама, Н. Метрополиса, Н. И. Бусленко, Дж. М. Хеммерсли, Дж. Спанье, И. М. Соболя, С. М. Ермакова, Г. А. Михайлова, Г. И. Марчука, М. Кейлоса и др. В последние шестьдесят лет область применения численного статистического моделирования значительно расширилась. Была разработана содержательная теория вероятностных представлений решений задач математической физики. На базе этой теории были построены эффективные(экономичные) оцениватели метода Монте-Карло. Алгоритмы численного статистического моделирования конструируются и используются также для решения задач статистической физики (схема Метрополиса — Хастингса, модель Изинга и др.), физической и химической кинетики (многочастичные схемы, численные методы решения уравнений Больцмана и Смолуховского, моделирование реакций и фазовых переходов, построение вероятностных клеточных автоматов и др.), теории массового обслуживания (марковские модели с пуассоновскими потоками заявок), финансовой математики (вероятностные формулы и итерационные процессы, алгоритмы решения стохастических дифференциальных уравнений и др.), теории турбулентности, математической биологии и т. д. Особое развитие получают также численные модели случайных процессов и полей (с применениями в метеорологии, физике атмосферы и океана и других областях), а также смешанные рандомизированные проекционные и дискретностохастические численные схемы (включая функциональные алгоритмы). Основная часть Общая схема метода Монте-Карло. Понятие монте-карловской оценки В самом общем виде схема метода Монте-Карло выглядит следующим образом. Пусть требуется приближенно вычислить на компьютере некоторую величину ${\displaystyle I}$. Предполагается, что имеется или строится (выбрается) случайная величина ${\displaystyle \zeta }$, математическое ожидание ${\displaystyle {\bf {E}}\zeta }$ которой равно ${\displaystyle I}$, дисперсия ${\displaystyle {\bf {D}}\zeta ={\bf {E}}(\zeta -{\bf {E}}\zeta )^{2}}$ конечна, и кроме того, выборочные значения ${\displaystyle \zeta _{i}}$ случайной величины ${\displaystyle \zeta }$ могут быть достаточно эффективно (экономично) численно смоделированы (реализованы на компьютере). АЛГОРИТМ. Численно моделируем (реализуем на компьютере) достаточно большое количество ${\displaystyle n}$  выборочных значений ${\displaystyle \zeta _{1},...,\zeta _{n}}$ случайной величины ${\displaystyle \zeta }$ и, используя закон больших чисел, получаем приближение требуемой величины: ${\displaystyle I={\bf {E}}\zeta \approx Z_{n}={\frac {S_{n}}{n}}={\frac {\zeta _{1}+...+\zeta _{n}}{n}}.}$ В русскоязычной литературе по методам Монте-Карло базовая случайная величина ${\displaystyle \zeta }$ называется оценкой величины ${\displaystyle I}$ (в англоязычной литературе используется термин «estimator», реже — «estimate»). Чаще всего монте-карловская оценка ${\displaystyle \zeta }$ имеет вид ${\displaystyle \zeta =q({\boldsymbol {\xi }}),}$ где ${\displaystyle {\boldsymbol {\xi }}=\left(\xi ^{(1)},...,\xi ^{(d)}\right)}$ — это случайный вектор или отрезок последовательности случайных величин с заданными плотностями распределениями ${\displaystyle f_{\boldsymbol {\xi }}({\bf {x}});\ {\bf {x}}\in X\subset \mathbb {R} ^{d},}$ а ${\displaystyle q({\bf {x}})}$ — функция (неслучайная) ${\displaystyle d}$ переменных ${\displaystyle {\bf {x}}=\left(x^{(1)},...,x^{(d)}\right)}$ (здесь ${\displaystyle d\leq \infty }$); при этом ${\displaystyle {\bf {E}}\zeta ={\bf {E}}q({\boldsymbol {\xi }})=\int _{X}q({\bf {x}})f_{\boldsymbol {\xi }}({\bf {x}})d{\bf {x}}.}$ Средняя погрешность и трудоемкость метода Монте-Карло В силу центральной предельной теоремы для достаточно большого числа смоделированных на компьютере выборочных значений ${\displaystyle \{\zeta _{1},...,\zeta _{n}\}}$ (которые можно трактовать как независимые одинаково распределенные случайные величины) имеем ${\displaystyle {\bf {E}}\delta _{n}={\bf {E}}|I-Z_{n}|={\bf {E}}\left|{\frac {S_{n}-n{\bf {E}}\zeta }{n}}\right|={\sqrt {\frac {{\bf {D}}\zeta }{n}}}\times {\bf {E}}\left|{\frac {S_{n}-n{\bf {E}}\zeta }{{\sqrt {{\bf {D}}\zeta }}{\sqrt {n}}}}\right|\approx }$ ${\displaystyle \approx {\bf {E}}\left|\xi ^{(0,1)}\right|\times {\sqrt {\frac {{\bf {D}}\zeta }{n}}}={\sqrt {\frac {2}{\pi }}}\times {\sqrt {\frac {{\bf {D}}\zeta }{n}}}\approx 0,80\times {\sqrt {\frac {{\bf {D}}\zeta }{n}}};}$ здесь ${\displaystyle \xi ^{(0,1)}}$ — стандартная гауссовская случайная величина с плотностью распределения ${\displaystyle f_{\xi ^{(0,1)}}(u)={\frac {e^{-u^{2}/2}}{\sqrt {2\pi }}};\;-\infty <u<+\infty }$. Таким образом, скорость сходимости метода Монте-Карло относительно невысока (это определяет основной недостаток метода). Чтобы получить следующий десятичный знак в мантиссе величины ${\displaystyle I}$ (т.е. уменьшить примерно в десять раз среднюю погрешность ${\displaystyle {\bf {E}}\delta _{n}}$), нужно численно смоделировать в сто раз больше выборочных значений ${\displaystyle \{\zeta _{i}\}}$. Поэтому в практических расчетах число ${\displaystyle n}$ берется весьма большим, а учитывая тот факт, что получение каждого ${\displaystyle \{\zeta _{1},...,\zeta _{n}\}}$ как правило, является трудоемким, можно констатировать, что методы Монте-Карло являются весьма «тяжеловесными» (затратными). Между тем важным преимуществом метода Монте-Карло является относительная простота подсчета затрат (а это, в свою очередь, следствие простоты основной схемы метода), что существенно упрощает применение естественного критерия оптимизации метода: тот алгоритм лучше, который позволяет достигнуть заданного уровня ${\displaystyle \Delta }$ средней погрешности за меньшее время. Вычислительные затраты метода Монте-Карло близки к времени вычисления суммы ${\displaystyle S_{n}}$ и приблизительно равны величине ${\displaystyle s=t\times n}$, где ${\displaystyle t}$ — это среднее время подсчета одного выборочного значения ${\displaystyle \zeta _{i}}$ случайной величины ${\displaystyle \zeta =q({\boldsymbol {\xi }})}$. Заметим, что ${\displaystyle \Delta }$ дисперсия ${\displaystyle {\bf {D}}\zeta }$ пропорциональна числу выборочных значений ${\displaystyle n}$. Действительно, из равенства ${\displaystyle 0,80\times {\sqrt {\frac {{\bf {D}}\zeta }{n}}}=\Delta \;\;}$следует соотношение ${\displaystyle \;\;n=\left({\frac {0,80}{\Delta }}\right)^{2}\times {\bf {D}}\zeta .}$ Это обосновывает возможность использования вместо затрат ${\displaystyle s}$ пропорциональной ей, более практичной и информативной величины ${\displaystyle S=t\times {\bf {D}}\zeta .}$ Эта величина называется трудоемкостью метода Монте-Карло и используется для формулировки следующего аналога сформулированного выше критерия качества вычислительного алгоритма: в алгоритмах метода Монте-Карло лучшим (оптимальным) является тот выбор случайной величины ${\displaystyle \zeta =q({\boldsymbol {\xi }})}$ и/или плотности распределения ${\displaystyle f_{\boldsymbol {\xi }}({\bf {x}})}$ вектора ${\displaystyle {\boldsymbol {\xi }}}$ , который дает меньшее значение трудоемкости ${\displaystyle S}$. Величины ${\displaystyle t}$ и ${\displaystyle {\bf {D}}\zeta }$ в выражении для трудоемкости оценивают с помощью предварительных расчетов для относительно малого числа выборочных знгачений ${\displaystyle {\hat {n}}\ll n}$. Методический пример: вычисление интеграла методом Монте-Карло Наиболее известным примером использования метода Монте-Карло является задача приближенного вычисления многократного интеграла ${\displaystyle I=\int _{X}g({\bf {x}})\,d{\bf {x}};\;\;\;X\subseteq \mathbb {R} ^{d}.}$ Выберем плотность распределения ${\displaystyle f_{\boldsymbol {\xi }}({\bf {x}}),\;{\bf {x}}\in \mathbb {R} ^{d},}$ случайного вектора ${\displaystyle {\boldsymbol {\xi }}}$ так, что ${\displaystyle f_{\boldsymbol {\xi }}({\bf {x}})>0}$ для ${\displaystyle g({\bf {x}})\neq 0}$ и выборочные значения ${\displaystyle {\boldsymbol {\xi }}_{1},\ldots ,{\boldsymbol {\xi }}_{n}}$ могут быть достаточно просто смоделированы на компьютере. Получить приближение интеграла ${\displaystyle I}$ можно следующим образом: ${\displaystyle I=\int _{X}g({\bf {x}})\,d{\bf {x}}=\int _{X}{\frac {g({\bf {x}})}{f_{\boldsymbol {\xi }}({\bf {x}})}}f_{\boldsymbol {\xi }}({\bf {x}})d{\bf {x}}={\bf {E}}\zeta ^{(I)}\approx Z_{n}^{(I)}={\frac {1}{n}}{\Biggl (}{\frac {g({\boldsymbol {\xi }}_{1})}{f_{\boldsymbol {\xi }}}}({\boldsymbol {\xi }}_{1})+\ldots +{\frac {g({\boldsymbol {\xi }}_{n})}{f_{\boldsymbol {\xi }}({\boldsymbol {\xi }}_{n})}}{\Biggr )},}$ где ${\displaystyle \zeta ^{(I)}=q({\boldsymbol {\xi }})=g({\boldsymbol {\xi }})/f_{\boldsymbol {\xi }}({\boldsymbol {\xi }})}$ весовая монте-карловская оценка интеграла ${\displaystyle I}$, а ${\displaystyle {\boldsymbol {\xi }}_{1},\ldots ,{\boldsymbol {\xi }}_{n}}$ — это выборочные значения случайного вектора ${\displaystyle {\boldsymbol {\xi }},}$ смоделированные на компьютере согласно выбранной плотности ${\displaystyle f_{\boldsymbol {\xi }}({\bf {x}}).}$ Средняя погрешность этого приближения равна ${\displaystyle {\bf {E}}\left|I-Z_{n}^{(I)}\right|\approx 0,80\times {\sqrt {\frac {{\bf {D}}\zeta ^{(I)}}{n}}},\;}$ где ${\displaystyle \;\;{\bf {D}}\zeta ^{(I)}=\int _{X}{\frac {g^{2}({\bf {x}})}{f_{\boldsymbol {\xi }}({\bf {x}})}}\,\,d{\bf {x}}-I^{2}.}$ Для сравнения заметим, например, что погрешность вычисления однократного интеграла ${\displaystyle I=\int _{a}^{b}g(x)\,dx}$ (здесь ${\displaystyle d=1}$) с гладкой подынтегральной функцией ${\displaystyle g(x)}$ с помощью простейшего метода центральных прямоугольников ${\displaystyle I\approx h\times \sum _{i=1}^{n}g\left({\frac {x_{i}+x_{i+1}}{2}}\right);\;\;x_{i}=a+(i-1)h;\;\;h={\frac {b-a}{n}}}$ (это вариант интегральной суммы) определяется соотношением ${\displaystyle n^{-2}}$ (на четыре порядка лучше, чем у алгоритма метода Монте-Карло). Погрешность чуть более «хитрой» формулы Симпсона Метод прямоугольников и формула Симпсона отражают частные случаи квадратурных формул Ньютона—Котеса, которые строятся с помощью соответствующих полиномиальных сеточных приближений подынтегральной функции ${\displaystyle g(x)}$. Хорошо известно, что при переходе к многомерному случаю и при рассмотрении негладких подынтегральных функций ${\displaystyle g({\bf {x}})}$ построение «хорошей» (в смысле скорости сходимости) сеточных кубатурных формул существенно затрудняется. В свою очередь, порядок сходимости ${\displaystyle n^{-1/2}}$ метода Монте-Карло не зависит от размерности ${\displaystyle d}$. Свойства подынтегральной функции ${\displaystyle g({\bf {x}})}$ (в том числе гладкость) влияют только на величину дисперсии ${\displaystyle {\bf {D}}\zeta }$. Поэтому «нишей» алгоритмов численного статистического моделирования считаются многомерные задачи со сложными (например, негладкими) входными данными (а таких задач на практике довольно много, если не большинство). В частности, при численном приближении интеграла ${\displaystyle I=\int _{X}g({\bf {x}})\,d{\bf {x}},\;X\subseteq \mathbb {R} ^{d}}$, считается, что — для размерностей ${\displaystyle 1\lesssim d\lesssim 3}$ целесообразно использовать детерминированные (сеточные) кубатурные формулы; — для размерностей ${\displaystyle d\gtrsim 10}$ (включая задачи бесконечной размерности) практически не имеет конкурентов метод Монте-Карло; —для размерностей ${\displaystyle 3\lesssim d\lesssim 10}$ можно использовать смешанные дискретно-стохастические методы численного интегрирования. Существенным преимуществом метода Монте-Карло для задачи вычисления интеграла является разнообразие возможностей уменьшать трудоемкость ${\displaystyle S=t\times {\bf {D}}\zeta }$, основной из которых является выборка по важности, основанная на том, что для двух вариантов выбора плотности — ${\displaystyle f_{{\boldsymbol {\xi }}^{(1)}}({\bf {x}})}$ и ${\displaystyle f_{{\boldsymbol {\xi }}^{(2)}}({\bf {x}})}$ — с примерно равными временами ${\displaystyle t}$ для реализации выборочных значений ${\displaystyle \left\{\zeta ^{(I,1)}=q^{(1)}({\boldsymbol {\xi }}^{(1)})=g({\boldsymbol {\xi }}^{(1)})/f_{{\boldsymbol {\xi }}^{(1)}}({\boldsymbol {\xi }}^{(1)})\right\}}$ и ${\displaystyle \left\{\zeta ^{(I,2)}=q^{(2)}({\boldsymbol {\xi }}^{(2)})=g({\boldsymbol {\xi }}^{(2)})/f_{{\boldsymbol {\xi }}^{(2)}}({\boldsymbol {\xi }}^{(2)})\right\}}$ тот вариант даст меньшую дисперсию (а значит, и трудоемкость ${\displaystyle S}$), для которого выбранная плотность ближе к функции ${\displaystyle H|g({\bf {x}})|}$ , где ${\displaystyle H=1{\Big /}\int _{X}|g({\bf {x}})|\,d{\bf {x}}}$. Компьютерное моделирование случайных величин, векторов и функций Из представленных выше материалов следует, что ключевым при реализации алгоритмов метода Монте-Карло является экономичное компьютерное моделирование (численная реализация выборочных значений ${\displaystyle {\boldsymbol {\xi }}_{i})}$ ${\displaystyle d}$-мерного случайного вектора ${\displaystyle {\boldsymbol {\xi }}=\left(\xi ^{(1)},...,\xi ^{(d)}\right)}$ с заданной совместной плотностью распределения ${\displaystyle f_{\boldsymbol {\xi }}({\bf {u}}).}$ Выбирается одно из ${\displaystyle d!}$ разложений функции ${\displaystyle f_{\boldsymbol {\xi }}({\bf {u}})}$ в произведение условных плотностей для компонента вектора (одномерных случайных величин) и происходит последовательное компьютерное моделирование этих компонент, причем на каждом шаге в условия условных плотностей подставляются выборочные значения, полученные на предыдущих шагах. В свою очередь, для моделирования (получения выборочных значений ${\displaystyle \xi _{i}}$) одномерных случайных величин ${\displaystyle \xi }$ выбирается один из подходящих (по точности и экономичности) алгоритмов, использующий стандартные случайные числа ${\displaystyle \alpha _{i}\in U(0,1)}$ (см. далее раздел 8). Здесь различают алгоритмы для дискретных случайных величин с распределением ${\displaystyle {\bf {P}}\{\xi =x_{i}\}=p_{i};\;p_{i}>0;\;\sum _{i=1}^{M}p_{i}=1;\;i=1,...,M;\;M\leq \infty }$ и непрерывных (точнее, абсолютно непрерывных) величин с заданной плотностью ${\displaystyle f_{\xi }(u)}$ на интервале ${\displaystyle (a,b);\;-\infty \leq a<b\leq +\infty }$ (здесь ${\displaystyle f_{\xi }(u)>0}$ при ${\displaystyle u\in (a,b);\;\int _{a}^{b}f_{\xi }(u)\,du=1}$ и ${\displaystyle f_{\xi }(u)=0}$ при ${\displaystyle u\notin (a,b))}$. Несложно построить аналоги (или комбинации) упомянутых ниже алгоритмов и для смешанных «дискретно-непрерывных» случайных величин и/или для непрерывных случайных величин, распределенных на нескольких отстоящих друг от друга интервалах, но эти случаи являются крайне редкими, «экзотическими». Для компьютерного моделирования дискретных и непрерывных одномерных случайных величин можно выделить: — стандартные алгоритмы и их модификации; — альтернативные алгоритмы широкого применения; — специальные алгоритмы. Стандартным алгоритмом моделирования дискретной случайной величины ${\displaystyle \xi }$ является такой, в котором принимается значение ${\displaystyle \xi _{0}=x_{m}}$, где номер ${\displaystyle m}$ получается вычитанием из реализованного на компьютере стандартного случайного числа ${\displaystyle \alpha _{0}\in U(0,1)}$ сумм вероятностей ${\displaystyle \sum _{j=1}^{m^{\prime }}p_{j}}$ до первого получения такого номера ${\displaystyle m^{\prime }=m}$, для которого ${\displaystyle \alpha _{0}-\sum _{j=1}^{m}p_{j}<0}$. Отметим также такую модификацию описанного стандартного алгоритма, как квантильный метод, весьма эффективный для случая, когда число ${\displaystyle M}$ значений дискретной случайной величины ${\displaystyle \xi }$ достаточно велико (вплоть до ${\displaystyle M=+\infty }$). Стандартным алгоритмом моделирования непрерывной случайной величины ${\displaystyle \xi }$ является метод обратной функции распределения, основанный на применении формулы ${\displaystyle \xi _{0}=F_{\xi }^{-1}(\alpha _{0});\;\;\alpha _{0}\in U(0,1);}$ здесь ${\displaystyle F_{\xi }(x)=\int _{-\infty }^{x}f_{\xi }(u)\,du}$ — функция распределения случайной величины ${\displaystyle \xi }$, монотонно возрастающая (а значит, имеющая обратную функцию) на интервале распределения ${\displaystyle (a,b)}$. При заданной плотности распределения ${\displaystyle f_{\xi }(u)}$ для получения вычислимой (т.е. представляющей собой композицию элементарных функций от стандартного случайного числа ${\displaystyle \alpha _{0}}$) формулы метода обратной функции распределения нужно разрешить уравнение ${\displaystyle \int _{a}^{\xi _{0}}f_{\xi }(u)\,du=\alpha _{0}}$ относительно верхнего предела интегрирования ${\displaystyle \xi _{0}}$ в элементарных функциях. Если такое решение удается получить, то плотность ${\displaystyle f_{\xi }(u)}$ называется элементарной. Примером элементарной плотности является плотность экспоненциального распределения ${\displaystyle f_{\xi ^{(\lambda )}}(u)=\lambda e^{-\lambda u};\;u>0,\lambda >0}$ (соответствующая формула обратной функции распределения: ${\displaystyle \xi _{0}^{(\lambda )}=-{\frac {\ln \alpha _{0}}{\lambda }})}$, а вот плотность стандартного гауссовского (нормального) распределения ${\displaystyle f_{\xi ^{(0,1)}}(u)={\frac {e^{-u^{2}/2}}{\sqrt {2\pi }}};\;-\infty <u<+\infty }$ не является элементарной. Литература Metropolis N., Ulam S. The Monte Carlo method // Journal of American Statistical Association. 1949. V. 44. No 249. P. 335–341. Бусленко Н. П., Голенко Д. И., Соболь И. М., Срагович В. Г., Шрейдер Ю. А. Метод статистических испытаний (метод Монте-Карло). М.: Физматгиз, 1962. Hammmersley J. M., Handscomb D. C. Monte Carlo Methods. New York: Jonh Wiley and Sons, 1964. Spanier J., Gelbard E. Monte Carlo Principles and Newtron Transport Problems. Addison–Wesley, Reading, 1969. Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1973. Михайлов Г. А. Некоторые вопросы теории методов Монте-Карло. Новосибирск: Наука, 1974. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1974. Марчук Г. И., Михайлов Г. А., Назаралиев М. А., Дарбинян Р. А., Каргин Б. А., Елепов Б. С. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М.: Наука, 1982. Kalos M. H., Whitlock P. A. Monte Carlo Methods. New York: Jonh Wiley & Sons, 1986. Михайлов Г. А. Оптимизация весовых методов Монте-Карло. М.: Наука, 1988. Сабельфельд К. К. Методы Монте-Карло в краевых задачах. М.: Наука, 1989. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло. М.: Академия, 2006. Войтишек А. В. Лекции по численным методам Монте-Карло. Новосибирск: НГУ, 2018. Михайлов Г. А., Войтишек А. В. Статистическое моделирование. Методы Монте-Карло. М.: Юрайт, 2025.