Метод Монте-Карло и его точность / forpes.ru

Главная
Метод Монте-Карло и его точность

Метод Монте-Карло и его точность +16

12.01.2016 17:21

vasiatka 9 11400 Источник

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

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

Метод имеет две основных особенности. Первая — простая структура вычислительного алгоритма. Вторая — ошибка вычислений, как правило, пропорциональна
$\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$ .

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

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

$P\left( \left| \frac{{{\rho }_{N}}}{N}-m \right|\le 3\frac{b}{\sqrt{N}} \right)\approx 0.9973.$
Если требуется конкретный уровень надежности вычислений $\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 )$ .

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

Нетрудно видеть, что линейная плотность лучше соответствует функции $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 с.

Комментарии (9)

RolexStrider
12.01.2016 21:33
#8736115
+2
Вот в одном абзаце у вас есть два по отдельности истинных утверждения:

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

В стандартных датчиках случайных чисел период чаще всего не превосходит 2 в степени разрядности операционной системы, а то и еще меньше.
… но из этого можно сделать ложный вывод, что «для вычислений по методу Монте-Карло чаще всего применяют стандартные датчики случайных чисел, период которых чаще всего не превосходит 2 в степени разрядности операционной системы, а то и еще меньше».
В действительности же применяются статистически-качественные ГПСЧ, чаще всего что-либо из семейства XorShift (период от 2^128-1 до 2^1024-1) либо Mersenne Twister — там уже от 2^19937 и более. А стандартные Math.rand(), которые в 85% случаев есть вариации на тему линейных конгруэнтных генераторов, да в каких-либо серьёзных симуляциях… Шутите? Вы бы еще RANDU вспомнили!
1. vasiatka
  12.01.2016 22:42
  #8736193
  Спасибо, за комментарий. Полностью с вами согласен. Последние замечания в статье направлены больше на студенческую аудиторию. Например, я в студенческие годы (лет 13 назад) был удивлен результатом моделирования (при N>10^10).
  Добавлю, что вариаций на тему генераторов случайных чисел очень много. Например, я часто использую для таких целей хорошо изученные криптографические алгоритмы (например, наши ГОСТы).
  1. RolexStrider
    12.01.2016 23:42
    #8736263
    +2
    Вспомнилось… В студенческие годы, и кстати примерно 13 лет назад напоролся на Modulo Bias (см.там соотв.параграф). Очень долго тупил: не мог понять, а что же я сделал не так с этим Math.rand() % m, ну типа всё очевидно же, получаем случайное число от 0...n, а нам надо 0...m, но ежели m < n, то поделим нацело, возьмем остаток, быстро, качественно, недорого, ~~херак-херак и в продакшн~~. Результат был, как бы немного предсказуем (см.статью по ссылке)… И неважно, взял бы я тогда стандартный сишный rand (он кстати и был), MT19937 или детектор космических лучей.
    Держу пари, что на этом «Modulo Bias» при симуляциях методом Монте-Карло гораздо больше народу споткнулось, чем на использовании негодных ГПСЧ.
    
    dimview
    13.01.2016 18:18
    #8737769
    +2
    Аналогичные грабли возникают, когда нужны случайные числа в плавающей точке, а ГПСЧ возвращает целые числа. Появляется желание просто поделить на RAND_MAX, но при этом есть неожиданные тёмные углы.

mickvav
13.01.2016 15:50
#8737517
+1
А чего MISER и VEGAS забыли?
Да, и оставлю тут
www.gnu.org/software/gsl/manual/html_node/Monte-Carlo-Integration.html
1. vasiatka
  13.01.2016 20:14
  #8737955
  Указанные алгоритмы несколько выходят за рамки данной заметки.
  Vegas — алгоритм или лучше сказать, конкретное решение, приближающее распределение используемой случайной величины к |g(x)| / I. Иными словами выбирает плотность распределения пропорционально |g(x)| на сколько это возможно. О необходимости этого я говорил.
  MISER — алгоритм, разбивает область интегрирования на две части рекурсивно до некоторой глубины. К каждой полученной области, применяется метод Монте-Карло. А количество итераций алгоритма в каждой области выбирается так, чтобы суммарная дисперсия получаемой оценки была минимальной.
  Но, возможно, ссылка на библиотеку кому-то пригодится.

dimview
13.01.2016 18:05
#8737745
+1
> и построить свой генератор

Зачем изобретать велосипед, когда есть PCG32? Ну или можно взять любой готовый блочный шифр и шифровать им числа от 0 до обеда.
1. vasiatka
  13.01.2016 21:30
  #8738035
  > и построить свой генератор
  
  Да, не удачно выразился. Не о велосипеде речь. Хотел сказать, что-то типа «найти подходящий генератор (проверенный, с нужными свойствами), отличный от стандартного, и реализовать его в своей программе», а получилось, как всегда.
  
  Зачем изобретать велосипед, когда есть PCG32?
  
  Я не хотел совсем касаться темы конкретных генераторов ПСЧ, хотелось только обратить внимание на важность их выбора.
  Что касается конкретного генератора.
  Тут выбор богатый, и PCG32 — это всего лишь один из многих. Если бы в решении этой задачи (генерации ПСЧ) все было бы так просто, то было бы ровно одно решение.
  Каждый выбирает, то что ему подходит лучше для решения конкретной задачи. Это мое мнение, возможно оно ошибочно, но пока такая концепция меня не подводила.

dimview
13.01.2016 21:52
#8738049
> Если бы в решении этой задачи (генерации ПСЧ) все было бы так просто,

Сейчас там действительно уже всё просто — если не известно, какой алгоритм использует стандартная библиотека, то надо на всякий случай взять готовую реализацию любого современного алгоритма (Mersenne Twister, PCG32) или блочного шифра (ChaCha20, Speck). Они все гарантированно проходят все тесты и отличаются только сложностью реализации и скоростью работы.

Преимущество PCG32 и Speck в том, что их реализация укладывается в несколько строчек, поэтому не нужно тянуть новую зависимость.

> то было бы ровно одно решение.

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