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

В самом методе нет ничего сложного. Именно эта простота объясняет популярность данного метода.

Метод имеет две основных особенности. Первая — простая структура вычислительного алгоритма. Вторая — ошибка вычислений, как правило, пропорциональна
\sqrt{D\zeta/N}, где D\zeta — некоторая постоянная, а N — число испытаний. Ясно, что добиться высокой точности на таком пути невозможно. Поэтому обычно говорят, что метод Монте-Карло особенно эффективен при решении тех задач, в которых результат нужен с небольшой точностью.

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



Общая схема метода


Допустим, что нам требуется вычислить какую-то неизвестную величину m. Попытаемся придумать такую случайную величину \xi, чтобы M\xi =m. Пусть при этом D\xi ={{b}^{2}}.
Рассмотрим N независимых случайных величин {{\xi }^{1}},{{\xi }^{2}},\ldots,{{\xi }^{N}} (реализаций), распределения которых совпадают с распределением \xi. Если N достаточно велико, то согласно центральной предельной теореме распределение суммы {{\rho }_{N}}=\sum\limits_{i}{{{\xi }_{i}}} будет приблизительно нормальным с параметрами M\rho_N=Nm, D\rho_N=Nb^2.

На основе Центральной предельной теоремы (или если хотите предельной теоремы Муавра-Лапласа) не трудно получить соотношение:

P\left( \left| \frac{{{\rho }_{N}}}{N}-m \right|\le k\frac{b}{\sqrt{N}} \right)=P\left( \left| \frac{1}{N}\sum\limits_{i}{{{\xi }_{i}}}-m \right|\le k\frac{b}{\sqrt{N}} \right)\to 2\Phi (k)-1,

где \Phi(x) — функция распределения стандартного нормального распределения.

Это — чрезвычайно важное для метода Монте-Карло соотношение. Оно дает и метод расчета m, и оценку погрешности.

В самом деле, найдем N значений случайной величины \xi. Из указанного соотношения видно, что среднее арифметическое этих значений будет приближенно равно m. С вероятностью близкой к (2\Phi (k)-1) ошибка такого приближения не превосходит величины kb/\sqrt{N}. Очевидно, эта ошибка стремится к нулю с ростом N.

В зависимости от целей последнее соотношение используется по разному:

  1. Если взять k=3, то получим так называемое «правило 3\sigma»:

    P\left( \left| \frac{{{\rho }_{N}}}{N}-m \right|\le 3\frac{b}{\sqrt{N}} \right)\approx 0.9973.

  2. Если требуется конкретный уровень надежности вычислений \alpha,
    P\left( \left| \frac{{{\rho }_{N}}}{N}-m \right|\le \left( {{\Phi }^{-1}}\left( (1+\alpha )/2 \right) \right)\frac{b}{\sqrt{N}} \right)\approx \alpha



Точность вычислений


Как видно из приведенных выше соотношений, точность вычислений зависит от параметра N и величины b – среднеквадратичного отклонения случайной величины \xi.

В этом пункте хотелось бы указать важность именно второго параметра b. Лучше всего это показать на примере. Рассмотрим вычисление определенного интеграла.

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

Итак, требуется вычислить определенный интеграл:

I=\int\limits_{a}^{b}{g(x)dx}

Выберем произвольную случайную величину \xi с плотностью распределения {{p}_{\xi }}(x), определенной на интервале (a,b). И рассмотрим случайную величину \zeta =g(\xi )/{{p}_{\xi }}(\xi ).

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

M\zeta =\int\limits_{a}^{b}{[g(x)/{{p}_{\xi }}(x)]{{p}_{\xi }}(x)dx=I}

Таким образом, получаем:

P\left( \left| \frac{1}{N}\sum\limits_{i}{{{\zeta }_{i}}}-I \right|\le 3\sqrt{\frac{D\zeta }{N}} \right)\approx 0.9973.

Последнее соотношение означает, что если выбрать N значений {{\xi }^{1}},{{\xi }^{2}},\ldots,{{\xi }^{N}}, то при достаточно большом N:

\frac{1}{N}\sum\limits_{i}{\frac{g({{\xi }_{i}})}{{{p}_{\xi }}({{\xi }_{i}})}\approx I}
.

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

Можно показать, что D\zeta будет иметь минимальное значение, когда {{p}_{\xi }}(x) пропорционально |g(x)|. Выбрать такое значение {{p}_{\xi }}(x) в общем случае очень сложно (сложность эквивалентна сложности решаемой задачи), но руководствоваться этим соображением стоит, т.е. выбирать распределение вероятностей по форме схожее с модулем интегрируемой функции.

Численный пример


Теория, конечно, дело хорошее, но давайте рассмотрим численный пример: a=0; b=\pi/2; g(x)=cos(x).

Вычислим значение интеграла с применением двух различных случайных величин.

В первом случае будем использовать равномерно распределенную случайную величину на [a,b], т.е. {{p}_{\xi }}(x)=2/\pi.

Во втором случае возьмем случайную величину с линейной плотностью на [a,b], т.е. {{p}_{\xi }}(x)=\frac{4}{\pi }(1-2x/\pi ).

Вот график, указанных функций

image

Нетрудно видеть, что линейная плотность лучше соответствует функции g(x).

Код программы модельного примера в математическом пакете Maple
restart;
with(Statistics):
with(plots):
#исходные функции
g:=x->cos(x):
a:=0:
b:=Pi/2:
N:=10000:
#плотности распределений
p1:=x->piecewise(x>=a and x<b,1/(b-a)):
p2:=x->piecewise(x>=a and x<b,4/Pi-8*x/Pi^2):
#графики
plot([g(x),p1(x),p2(x)],x=a..b, legend=[g,p1,p2]);
#Точное значение интеграла
I_ab:=int(g(x),x=0..b);
#функция метода Монте-Карло для вычисления приближенного вычисления интеграла
#не стоит ее использовать при реальных расчетах
INT:=proc(g,p,N)
  local xi;
  xi:=Sample(RandomVariable(Distribution(PDF = p)),N);
  evalf(add(g(xi[i])/p(xi[i]),i=1..N)/N);
end proc:
#Приближенное значение интеграла
I_p1:=INT(g,p1,N);#c равномерной плотностью
I_p2:=INT(g,p2,N);#c линейной плотностью
#Абсолютная погрешность
Delta1:=abs(I_p1-I_ab);#c равномерной плотностью
Delta2:=abs(I_p2-I_ab);#c линейной плотностью
#Относительные погрешности в процентах
delta1:=Delta1/I_ab*100;#c равномерной плотностью
delta2:=Delta2/I_ab*100;#c линейной плотностью
#Вычисление дисперсий
Dzeta1:=evalf(int(g(x)^2/p1(x),x=a..b)-1);
Dzeta2:=evalf(int(g(x)^2/p2(x),x=a..b)-1);
#Оценка погрешности в первом случае
3*sqrt(Dzeta1)/sqrt(N);
#Оценка погрешности во втором случае
3*sqrt(Dzeta2)/sqrt(N);

Файл с данной программой можно взять тут


Точное значение интеграла легко вычислить аналитически, оно равно 1.

Результаты одного моделирования при N=10:

Для равномерно распределенной случайной величины: I\approx 1.21666.

Для случайной величины с линейной плотностью распределения: I\approx 0.97641.

В первом случае относительная погрешность более 21%, а во втором 2.35%. Точность 3\sqrt{\frac{D\zeta }{N}} в первом случае равна 0.459, а во втором – 0.123.

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

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

Количество итераций и генераторы случайных чисел


Не трудно видеть, что точность вычислений зависит от количества N случайных величин включенных в сумму. Причем, для увеличения точности вычислений в 10 раз нужно увеличить N в 100 раз.

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

Если в качестве источника случайности используется некоторое физическое явление (физический датчик случайных чисел), то все работает отлично.

Часто для вычислений по методу Монте-Карло применяют датчики псевдослучайных чисел. Главная особенность таких генераторов – наличие некоторого периода.

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

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

Литература


Популярные лекции по математике 1968. Выпуск 46. Соболь И.М. Метод Монте-Карло. М.: Наука, 1968. — 64 с.