Метод Монте-Карло и его точность / Хабр

Введение

Метод Монте-Карло – это численный метод решения математических задач при помощи моделирования случайных величин.

Датой рождение метода Монте-Карло принято считать 1949 г., когда появилась статья под названием «Метод Монте-Карло» (Н. Метрополис, С. Улам). Создателями этого метода считают американских математиков Дж. Неймана и С. Улама. В нашей стране первые статьи были опубликованы в 1955–56 гг. (В.В. Чавчанидзе, Ю.А. Шрейдер, В.С. Владимиров)

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

Само название «Монте-Карло» происходит от города Монте-Карло в княжестве Монако, знаменитого своим игорным домом, а одним из простейших механических приборов для получения случайных величин является рулетка.

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

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

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

1.Теоретическая часть

1.1 Сущность метода Монте-Карло и моделирование случайных величин

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

Рис. 1.1

Предположим, что эта фигура расположена внутри единичного квадрата.

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

Для того чтобы выбирать точки случайно, необходимо перейти к понятию случайная величина. Случайная величина непрерывная, если она может принимать любое значение из некоторого интервала .

Непрерывная случайная величина определяется заданием интервала , содержащего возможные значения этой величины, и функции , которая называется плотностью вероятностей случайной величины (плотностью распределения ). Физический смысл следующий: пусть – произвольный интервал, такой что , тогда вероятность того, что окажется в интервале , равна интегралу

(1.1)

Множество значений может быть любым интервалом (возможен случай ). Однако плотность должна удовлетворять двум условиям:

1) плотность положительна:

; (1.2)

2) интеграл от плотности по всему интервалу равен 1:

(1.3)

Математическим ожиданием непрерывной случайной величины называется число

(1.4)

Дисперсией непрерывной случайной величины называется число:

Нормальной случайной величиной называется случайная величина , определённая на всей оси и имеющая плотность

(1.5)

где – числовые параметры

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

, называемой обычно интегралом вероятностей.

Согласно (1.1)

В интеграле сделаем замену переменной , тогда получим

,

где Отсюда следует, что Также

Нормальные случайные величины очень часто встречаются при исследовании самых различных по своей природе вопросов.

Выбрав , , найдём . Следовательно,

(1.6)

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

Проводя большое количество опытов, и получая большое количество случайных величин можно воспользоваться центральной предельной теоремой теории вероятностей. Эта теорема впервые была сформулирована П. Лапласом. Обобщением этой теоремы занимались многие выдающиеся математики, в том числе П.Л. Чебышёв, А.А. Марков, А.М. Ляпунов. Её доказательство достаточно сложно.

Рассмотрим одинаковых независимых случайных величин , так что распределения вероятностей этих величин совпадают. Следовательно, их математические ожидания и дисперсии также совпадают. Величины эти могут быть как непрерывными, так и дискретными.

Обозначим

Сумму всех этих величин обозначим через

Используя соотношения

получаем

Рассмотрим теперь нормальную случайную величину с такими же параметрами: .

В центральной предельной теореме утверждается, что для любого интервала при больших

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

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

Рассмотрим независимых случайных величин распределения которых совпадают с распределением . Если достаточно велико, то, согласно центральной предельной теореме, распределение суммы будет приблизительно нормальным с параметрами . Из (1.6) следует, что .

Последнее соотношение перепишем в виде:

(1.7)

Это соотношение даёт и метод расчёта , и оценку погрешности.

В самом деле, найдём значений случайной величины . Из (1.7) видно, что среднеарифметическое этих значений будет приближенно равно . С большой вероятностью погрешность приближения не превосходит величины . Эта погрешность стремится к нулю с ростом . На практике часто используют не оценку сверху , а на вероятную ошибку, которая приближенно равна Именно такой обычно порядок фактической погрешности расчёта, которая равна

.

Для получения случайных чисел используют обычно три способа: таблицы случайных величин, генераторы случайных чисел и метод псевдослучайных чисел.

Таблицы случайных чисел используют предпочтительно при расчётах вручную. Определяющую роль в этом играют два факта: 1) при использовании ЭВМ легче и удобней воспользоваться генератором случайных чисел, получаемых тут же, чем загружать из памяти значения таблицы, которая к тому же, будет занимать там место. 2) При подсчёте вручную нет необходимости использовать ЭВМ, так как часто необходимо выяснить лишь порядок искомой величины.

Генераторы случайных чисел анализируют какой-либо процесс, доступный для них (шумы в электронных лампах, скачки напряжения) и составляют последовательность из 0 и 1, из которых составляются числа с определёнными разрядами, однако такой метод получения случайных величин имеет свои недостатки. Во-первых, трудно проверить вырабатываемые числа. Проверки приходится делать периодически, так как из-за каких-либо неисправностей может возникнуть так называемый дрейф распределения (нули и единицы в каком-либо из разрядов станут появляться не одинаково часто). Во-вторых, обычно все расчёты на ЭВМ проводятся несколько раз, чтобы исключить возможность сбоя. Но воспроизвести те же самые случайные числа невозможно, если их только не запоминать по ходу счёта. А если запоминать, то снова появляется случай таблиц.

Таким образом, самым эффективным способом получения случайных чисел – это использование псевдослучайных чисел.

Числа, получаемые по какой-либо формуле и имитирующие значения случайной величины , называются псевдослучайными числами.

Первый алгоритм для получения псевдослучайных чисел был предложен Дж. Нейманом. Он называется методом середины квадратов.

Пусть задано 4-значное число . Возведём его квадрат. Получим 8-значное число . Выберем 4 средние цифры этого числа и положим .Далее и т.д.

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

при (1.8)

По числам вычисляются псевдослучайные числа

(1.9)

Формула (1.8) означает, что число равно остатку, полученному при делении на , такой остаток называют наименьшим положительным вычетом по модулю Формулы (1.8), (1.9) легко реализовать на ЭВМ.

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

Единственный недостаток метода – ограниченность количества псевдослучайных чисел, так как если последовательность чисел вычисляется на ЭВМ по формуле вида

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

Значения любой случайной величины можно получить путём преобразования значений одной какой-либо случайной величины. Обычно роль такой случайной величины играет случайная величина , равномерно распределённая в . Процесс нахождения значения какой-либо случайной величины путём преобразования одного или нескольких значений называется разыгрыванием случайной величины .

Допустим, что необходимо получать значения случайной величины , распределённой в интервале , с плотностью . Докажем, что значения можно находить из уравнения

(1.10)

т.е. выбрав очередное значение , надо решить уравнение (1.10) и найти очередное значение .

Для доказательства рассмотрим функцию

.

Из общих свойств плотности (1.2), (1.3) следует, что

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

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

Поэтому, если принадлежит интервалу , то принадлежит интервалу , и наоборот. Значит

Так как равномерно распределена в , то

,

итак

,

а это и означает, что случайная величина , являющаяся корнем уравнения (1.10), имеет плотность вероятностей .

Может оказаться, что разрешить уравнение (1.10) относительно трудно, например, в случаях, когда интеграл от не выражается через элементарные функции или когда плотность задана графически. Предположим, что случайная величина определена на конечном интервале и плотность её ограничена .

Разыгрывать значение можно следующим образом:

1) выбираются два значения и случайной величины и строится случайная точка с координатами

2) если точка лежит под кривой , то полагаем , если же точка лежит над кривой , то пара отбрасывается и выбирается новое значение.

1.2 Вычисление интегралов

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

(2.1)

Этот интеграл может быть несобственным, но абсолютно сходящимся.

Выберем произвольную плотность распределения , определённую на интервале . Наряду со случайной величиной , определённой в интервале с плотностью , необходимо определить случайную величину

Согласно соотношению получим

Рассмотрим теперь одинаковых независимых случайных величин и применим к их сумме центральную предельную теорему. Формула (1.7) в этом случае запишется так:

Последнее соотношение означает, что если выбирать значений , то при достаточно большом

(2.2)

Оно показывает также, что с очень большой вероятностью погрешность приближения (2.2) не превосходит .

Для расчёта интеграла (2.1) можно использовать любую случайную величину . Определённую в интервале с плотностью . В любом случае . Однако дисперсия , а с ней и оценка погрешности формулы (2.2) зависят от того, какая величина используется, так как

(2.3)

Докажем, что это выражение будет минимальным тогда, когда пропорциональна .

Для этого воспользуемся неравенством

, в которым положим , . Получим неравенство

(2.4)

Из (2.3), (2.4) следует, что

(2.5)

Остается доказать, что нижняя граница дисперсии (2.5) реализуется при выборе плотности . Так как

.

Следовательно,

,

и правая часть (2.3) обращается в правую часть (2.5)

Использовать плотность для расчёта практически невозможно, так как для этого нужно знать значение интеграла . А его вычисление представляет собой задачу, равноценную задаче о вычислении интеграла (2.1). Поэтому ограничиваются следующей рекомендацией: желательно, чтобы плотность была пропорциональна .

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

Также если стоит задача вычислить интеграл (2.1), преобразуем его к виду

(2.6)

Если теперь обозначить (2.7)

То интеграл принимает вид

(2.8)

и может быть вычислен при помощи метода статистических испытаний.

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

Как известно, плотность вероятности равномерного закона распределения в интервале равна:

(2.9)

Подставим в интеграл (2.6) значение из формулы (2.9) и получим:

(2.10)

и рассмотрим процедуру вычисления:

из множества равномерно распределённых случайных чисел выбирается . Для каждого значения вычисляется , затем вычисляется среднее значение

(2.11)

функции на интервале

Таким образом, величина интеграла (2.10) может быть представлена в виде следующей формулы

(2.12)

Рассмотренный частный случай находит широкое применение интегралов методом статистического моделирования в силу того, что границы области определения могут быть легко приведены к пределам интегрирования

1.3 Вычисление кратных интегралов

Обычно при вычислении кратных интегралов методом Монте-Карло используют один из двух способов.

Первый способ.

Пусть требуется вычислить кратный интеграл

(3.1)

по области G, лежащей в мерном единичном кубе

Выберем равномерно распределённых на отрезке последовательностей случайных чисел

Тогда точки можно рассматривать как случайные, равномерно распределённые в мерном единичном кубе.

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

(3.2)

где под понимается мерный объём области интегрирования. Если вычисление объёма затруднительно, то можно принять , и для приближенного вычисления интеграла получим:

(3.3)

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

Второй способ.

Если функция , то интеграл (3.1) можно рассматривать как объём тела в мерном пространстве, т.е.

(3.5)

где область интегрирования определяется условиями

Если в области , то введя новую переменную , получим

где область лежит в единичном мерном кубе

Возьмём равномерно распределенных на отрезке случайных последовательностей

Составим соответствующую последовательность случайных точек

Пусть из общего числа случайных точек точек принадлежат объёму , тогда имеет место приближенная формула

(3.6)

2. Практическая часть

2.1 Пример 1

Вычислим приближенно интеграл

Точное значение его известно:

Используем для вычисления две различные случайные величины , с постоянной плотностью (т.е. равномерна распределена в интервале ) и с линейной плотностью .Линейная плотность более соответствует рекомендации о пропорциональности и . Поэтому следует ожидать, что второй способ вычисления даст лучший результат.

1) Пусть , формула для разыгрывания имеет вид . А формула (2.2) примет вид .

Пусть . В качестве значений используем тройки чисел из табл. 1 (см. приложение), умноженные на 0.001. Промежуточные результаты сведены в табл. 2.1. Результат расчёта

Таблица 2.1

2) пусть теперь . Для разыгрывания используем формулу

,

откуда получаем

формула (2.2) имеет вид

Пусть . Числа выберем те же, что и в случае 1. Промежуточные результаты сведены в табл. 2.2. Результат расчёта

Таблица 2.2

Как и ожидалось, второй способ вычислений дал более точный результат.

3) По значениям, приведённым в таблицах (2.1) и (2.2) можно приближенно сосчитать дисперсии для обоих методов расчёта:

для 1:

для 2:

Несмотря на то, что значение невелико и приближенная нормальность оценки (2.2) не гарантирована, вычислим для обоих методов величины . Получим значения 0.103 и 0.027. Также фактические абсолютные погрешности при расчёте , равные 0.048 и 0.016, – величины того же порядка. Точные же значения в рассмотренном примере равны 0.233 и 0.0166. Таким образом, и при оценке дисперсий метод 2 оказался точнее метода 1.

2.2 Пример 2

Рассмотрим пример:

Требуется вычислить интеграл

(3.4)

где область G задаётся следующими неравенствами:

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

Записываем координаты и случайных точек в табл. 3.1, округляя до 3 знаков после запятой, и выбираем те из них, которые принадлежат области интегрирования.

Заполним табл. 3.1 по правилу:

1) Среди всех значений выделяем те, которые заключены между и .Для этих значений полагаем , для всех остальных

2) Среди всех значений . Соответствующих выделенным , выбираем те, которые заключены между

Для этих значений полагаем , для всех остальных

Таблица 3.1

3) Вычисляем . Области тнтегрирования принадлежат только те точки, для которых . В примере

4) Вычисляем значения подынтегральной функции в полученных точках.

После заполнения табл. 3.1 вычисляем площадь области интегрирования и по формуле (3.2) находим

Для сравнения приведём точное значение интеграла

Результат имеет сравнительно небольшую точность потому, что число точек недостаточно велико.

2.3 Пример 3

Рассмотрим пример: найдём приближенно объём, ограниченный поверхностями

Искомый объём численно равен величине интеграла

(3.7)

Так как в области V, вводим новую переменную , в результате чего интеграл (3.7) переходит в интеграл

(3.8)

где область, ограниченная поверхностями

т.е. принадлежит единичному кубу .

Берём теперь три равномерно распределенные на отрезке последовательности случайных чисел и записываем их в качестве координат случайных точек в табл. 3.2. Затем проверяем, какие из этих точек принадлежат области .

Таблица 3.2

Заполним табл. 3.2 по правилу:

1) выделяем точки, у которых , и полагаем для них

2) среди выделенных точек области принадлежат те, для которых выполняется неравенство .

Для этих точек , для остальных

3) вычисляем . Области принадлежат те точки, для которых

4) среди точек, у которых , области принадлежат те точки, координаты которых удовлетворяют неравенству

Для этих точек .

В примере общее количество точек , а число точек, принадлежащих области , равно 15. По формуле (3.6) получаем

, а точное значение объёма равно

Погрешность формулы (3.6) обратно пропорциональна корню из числа испытаний, т.е. .

Это означает, что для обеспечения большой точности число точек должно быть очень велико. Но так как приближенные формулы (3.3), (3.6) не зависят от размерности интеграла, метод Монте-Карло оказывается выгодным при вычислении интегралов большой размерности.

Заключение

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

Были рассмотрены основные свойства метода Монте-Карло и создана программа, показывающая возможности данного метода при использовании ЭВМ.

Было выяснено, что методом Монте-Карло можно решать разнообразные задачи, в том числе вычисление интегралов, не прибегая к сложным математическим вычислениям. Простота алгоритма метода Монте-Карло позволяет успешно реализовывать их на ЭВМ.

Список литературы

1. Бусленко Н.П. Метод статистического моделирования – М.: Статистика, 1970. – 112 с.

2. Демидович Б.П., Марон И.А. Основы вычислительной математики. – М.: Наука, 1966. – 664 с.

3. Епанешников А.М., Епанешников В.А. Программирование в среде TURBOPASCAL 7.0 – М.: Диалог-МИФИ, 1998. – 288 с.

4. Ермаков С.М. Метод Монте-Карло и смежные вопросы – М.: Наука, 1975–472 с.

5. Копченова Н.В., Марон И.А. Вычислительная математика в примерах и задачах. – М.: Наука, 1972. – 367 с.

6. Соболь И.М. Метод Монте-Карло – М.: Наука, 1985. – 80 c.

Приложения

1. Таблица 400 случайных цифр

2. Таблица 40 случайных чисел, равномерно распределенных на отрезке

3. Листинг программы

Вычисляются значения кратных интегралов из примера 2–3.

program pmk;

uses crt;

var

w, u, h, k, v, y, p, s, g, x, x2, y2, z2, niu, Integral, Integral2:real;

n, m, i, a, b, e1, e2, e, e3, e4, e5:integer;

begin

clrscr;

writeln (‘vychisleniye dvoynogo integrala iz primera 1’);

writeln (‘vvedite kolichestvo sluchaynykh tochek:’);

readln(n);

for i:=1 to n do

begin

g:=random;

p:=random;

x:=g;

y:=p;

if ((0.5<=x) and (x<=1)) then e1:=1

else e1:=0;

if ((0<=y) and (y<=2*x-1)) then e2:=1

else e2:=0;

e:=e1*e2;

if e=1 then s:=s x*x y*y;

if e=1 then a:=a 1;

v:=1/4;

delay(1000);

end;

Integral:=(v/a)*(s);

writeln (‘summa=’, s:5:5);

writeln (‘dvoynoy integral iz 1 primera =’, Integral:5:5);

writeln (‘vychisleniye troynogo integrala iz primera 2’);

writeln (‘vvedite kolichestvo sluchaynykh tochek:’);

readln(m);

for i:=1 to m do

begin

w:=random;

u:=random;

h:=random;

x2:=w;

y2:=u;

niu:=h;

if niu<=0.8 then e3:=1;

if (x2–0.5)*(x2–0.5) (y2–0.5)*(y2–0.5)<=(0.5)*(0.5) then e4:=1

else e4:=0;

e5:=e3*e4;

if (((0.8<niu) and (niu<1)) and ((x2–0.5)*(x2–0.5) (y2–0.5)*(y2–0.5) 6.25*(niu-0.8)*(niu-0.8)<=(0.5)*(0.5))) then e5:=1;

if e5=1 then b:=b 1;

delay(1000);

end;

Integral2:=2.5*(b/m);

writeln (‘kvo pod t =’, b:5);

writeln (‘troynoy integral iz 2 primera =’, Integral2:5:5);

readln;

end.

4. Пример работы программы при 10000 случайных точек