Раздел. 
Общие сведения про метод Монте-Карло

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

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

Этот метод позволяет решать задачи, в условиях которых присутствует элемент неопределенности (например, при подбрасывании монеты может выпасть «орел» или «решка»). Пусть требуется найти некоторую величину (например, долю выпадения «орлов»). На ЭВМ с помощью генератора (датчика) случайных чисел (ДСЧ) имитируются ситуации или процессы, возможные по условию задачи, и которые приводят к тем или иным исходам, при этом искомая величина принимает некоторое значение (в нашем примере это 0 — если выпал «орел» и 1 — в противном случае). Все или почти все различные исходы (с учетом, когда монета может упасть на ребро) проявятся, если многократно рассмотреть случайное развитие одного и того же начального состояния (смоделировать некоторое количество историй — N).

Закон больших чисел «разыгрываемых» историй утверждает, что среднее арифметическое полученных в каждом розыгрыше значений исследуемой величины имеет предельное искомое значение (в нашей задаче оно равно ½).

Это вероятностная сходимость, т. е. чем больше историй, тем достоверней можно утверждать, что наш результат близок к истинному. Для задач с элементами неопределенности — а в реальном мире все задачи такие — это даже естественно. Погрешность определения предельного значения пропорциональна 1/vN. Таким образом, для увеличения точности результата на один порядок, требуется разыграть в 100 раз больше историй. «Разумно» число историй первоначально задавать порядка 10 000.

Следовательно, требуется получать много случайных чисел так, чтобы переход от одного числа к другому определялся простыми правилами, но чтобы сами числа «производили впечатление случайности» (их называют поэтому псевдослучайными числами). Например, для выбора последовательности случайных цифр можно взять дробную часть числа пи (р = 3.14 159 265 358 979 322 710 605 500 359 114 752 4 197 169 399 375 106 094 117 768 666 367 115 791 941 615 428 947 141 787 123 712… Представляет интерес обратное утверждение: возможно ли, что для любой конечной, наперед заданной последовательности цифр есть ее вложение в бесконечном представлении числа р). Данный способ, однако, мало приемлем для программирования. Как правило, при решении задач методом Монте-Карло используются процедуры, которые с помощью рекуррентных формул генерируют случайные числа, равномерно распределенные на промежутке [0, 1]. В Pascal для этого используется стандартная функция RANDOM. Для отладки программы бывает важно уметь воспроизвести псевдослучайные числа, а для генерации другой последовательности случайных чисел используется процедура RANDOMIZE.

Очень полезным для понимания данного метода является моделирование игровых вероятностных ситуаций (бросание монеты или кубика, блуждания). Именно «игровая» или связанная с чем-то знакомым (известным, «бытовым») формулировка задачи помогает лучше усвоить метод, осмыслить понятие вероятности.

Пример 1. (Районная олимпиада).

Три игрока (с номерами 1, 2 и 3), имеющие изначально X, Y и Z жетонов соответственно, играют в следующую игру. В каждом раунде каждый игрок ставит на кон один жетон. Затем бросают кубик, на котором цифры 4, 5, 6 заменены на 1, 2 и 3. При выпадении числа i игрок с номером i забирает с кона все три жетона. Игра заканчивается, когда кто-нибудь из игроков проигрывает все жетоны. Введем функцию f (X, Y, Z), как среднюю длительность игры (среднее количество раундов) при заданных начальных капиталах X, Y, Z. Например, f (2, 2, 2) = 2. Ваша задача состоит в том, чтобы определить эту функцию. Для этого необходимо смоделировать игру на компьютере, накопить экспериментальные результаты, проанализировать их, а затем выдвигать гипотезы о виде функции f, проверять их для разных входных значений, и, отбросив неподходящие, найти решение.

Моделирование игры не вызывает трудности (программа приведена ниже). Также очевидно, что вид функции симметричен относительно порядка задания входных параметров f (X, Y, Z) = f (Y, X, Z) и т. д. Сложность задачи заключается в нахождении вида функции f (X, Y, Z) = X · Y · Z / (X + Y + Z — 2), так как результаты моделирования определяются не точно.

var.

i, n, Rounds: LongInt;

x, y, z, j: Integer;

a: array[0.2] of Integer;

begin.

Write ('x, y, z натуральные: '); ReadLn (x, y, z);

Write ('Число историй (игр): '); ReadLn (n);

Rounds:=0; {число раундов в одной игре}.

Randomize;

for i:=1 to n do {разыгрываем n историй развития начальной ситуации}.

begin.

a[0]: =x; a[1]: =y; a[2]: =z;

while a[0]*a[1]*a[2]0 do {моделируем игру, пока у всех игроков есть жетоны}.

begin.

Inc (Rounds);

for j:=0 to 2 do Dec (a[j]); {все игроки ставят по одному жетону}.

Inc (a[Random (3)], 3); {а победитель берет все три жетона}.

end;

end;

WriteLn ('f (x, y, z) = ', Rounds/n:0:3); {среднее количество раундов}.

end.

Данная задача достаточно хорошо характеризует метод Монте-Карло, а именно:

Идею метода.

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

То есть результат приближенно равен числу.

.

где x_i — результат игры i, а Nчисло всех проведенных игр (испытаний) Достоинство метода.

— незнание a priori (до опыта) функциональных зависимостей исследуемой задачи в целом, выявление этих зависимостей a posteriori (после опыта).

Недостатки метода.

• — неопределенное время расчета (варианты примера 1 при больших числах X, Y, Z);
• — приближенное вычисление результата.

Последний недостаток компенсируется тем, что с использованием данного метода вместе со значением может одновременно определяться и его погрешность по формулам:

При больших N формулу можно упростить:

где.

В пределах [ -, + ] с достоверностью 68.27% находится искомая величина, а в пределах ± 2 или ± 3 достоверность уже 95.45% и 99.73% соответственно. Поэтому метод по праву называют порой прецизионным или точным в смысле, что известна точность рассчитываемых величин, и это может служить точкой отсчета для проверки программ, использующих другие приближенные методы.

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

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

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

Задача 1. «Лучшее пари для простаков». (Районная олимпиада 1997).

Игрок A выбирает комбинацию из цифр 0 и 1 длиной 3 знака (например, 001). Игрок B выбирает свою комбинацию (отличную от игрока A). Подбрасывается монета и записываются результаты бросания (например, 101 101…, где 0 обозначает «орел», а 1 — «решка»). Игра прекращается в тот момент, когда в последовательности цифр на конце возникает комбинация, выбранная A или B (побеждает A или B соответственно). Игра повторяется.

• а) Оценить шансы на выигрыш каждого из игроков R (A, B) (т.е. отношение числа выигрышей игрока B к числу выигрышей игрока A).
• б) Для выбранной игроком A комбинации определить такую комбинацию для игрока B, которая ему дает больше шансов на выигрыш.

Ниже представлена таблица значений R (A, B) для всевозможных выбранных игроками A и B исходных комбинаций при «неограниченном продолжении» игры (выделены наиболее выигрышные ситуации для игрока B).

Пари является беспроигрышным (!) для игрока B. Парадокс заключается в том, что какую бы комбинацию цифр не выбрал игрок A, его соперник B может выбрать другую комбинацию, которая ему дает больше шансов на выигрыш.

Таблица R (A, B).

Указание. При решении задачи полученные результаты по пункту a) не будут совпадать с данными из таблицы, так как число опытов ограничено, тем не менее, позволяют дать качественный ответ по пункту б) Метод Монте-Карло используется для определения вероятности наступления какого-либо события.

Задача 2.

Пусть дана ось с отмеченными на ней целочисленными точками. Предположим, что Чиполлино спрятался в точке 0, в точке N находится пропасть, и сыщик Моркоу находится в точке k (0 < k < N). Сыщик ищет Чиполлино случайным образом, блуждая по соседним целочисленным точкам. Если он попадет в точку 0, то найдет Чиполлино, а если попадет в точку N, то свалится в пропасть. С какой вероятностью сыщик найдет Чиполлино?

Под вероятностью какого-либо события (P) мы будем понимать предельное значение частоты события, а именно, отношения числа успешных (приведших к появлению данного события) испытаний (N_у) к общему числу проведенных испытаний (N), то есть P «N_у / N. Чем больше мы проведем испытаний, тем точнее (в идеале) мы определяем численное значение вероятности. Очевидно, что вероятность P удовлетворяет условию: 0 Ј P Ј 1.

Указание. Смоделируйте многократный поиск сыщика из точки k, доля удачных попыток от общего их числа дает приближенную оценку искомой вероятности. На основании этой оценки сформулируйте простую формулу для нахождения вероятности события (обозначим ее P₀(k)), указанного в задаче.

Данную задачу можно решить точно, используя рекуррентную формулу:

P₀(k) = P₀(k-1)/2 + P₀(k+1)/2 при очевидных условиях, что P₀(0) = 1 и P₀(N)=0, однако использовать рекурсию в такой форме неприемлемо, так как глубина рекурсии неограниченна, что ведет к переполнению стека и краху программы.

Упражнение 1.

Выведите рекуррентную зависимость P₀(k+1) от P₀(k), начиная с k = 0 и, используя известное значение P₀(N), обратным ходом получите общее выражение для P₀(k).

Задача 3.

Пусть имеется «однорукий бандит» — игровой автомат с ручкой, которой его запускают для игры. Считаем (для упрощения), что игра будет типа «в орлянку», и игрок имеет начальный капитал в одну монету. Игра ведется до тех пор, пока игрок не обанкротится или выиграет Nмонет. Промоделировать игру для N=10, определить вероятность (шансы) игрока «сорвать указанный куш» и объяснить, почему автомат назвали «бандитом».

Задача о разорении игрока аналогична задачи блуждания по отрезку [0, N] с той лишь разницей, что требуется определить вероятность P_N(k) достичь точки N, находясь в точке k=1. Образно говоря, эти модели «связаны одной цепью». Последовательность испытаний, в которой каждое следующее испытание зависит только от исхода предыдущего, называется цепью Маркова. Многие реальные явления (например, броуновское движение частиц, обслуживание телефонных линий) описываются данными вероятностными моделями /2/.

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

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

Пример 2. (Областная олимпиада 2001).

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

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

Переформулируем условие задачи. Опишем квадрат около одной из окружностей (например, меньшего радиуса). Будем случайным образом кидать точки в этот квадрат. При достаточно большом их количестве ($n$) они равномерно распределятся по площади квадрата. Часть из них (m) попадет в область пересечения трех окружностей. Можно ожидать, что отношение m/$n$ имеет конечный предел, равный отношению искомой площади к площади описанного квадрата (см. Упражнение 2).

Указание. Наилучший путь — это «использовать геометрию» для анализа частных случаев (когда нет пересечения, одна окружность внутри другой), а метод МонтеКарло — для общего случая.

label.

NotInCircle;

var.

i, n, m: LongInt;

j: Integer;

x, y, r, rr: array[1.3] of Real;

xp, yp, xmin, ymin, d: Real;

begin.

for j:=1 to 3 do {вводим координаты центра и радиус трех окружностей}.

begin.

Write (j, '-я окружность (x y r): ');

ReadLn (x[j], y[j], r[j]);

rr[j]: =Sqr (r[j]); {вычислим квадраты радиусов — они будут часто использоваться}.

end;

Write ('Количество историй: '); ReadLn (n);

xmin:=x[1]- r[1]; {опишем квадрат около первой окружности}.

ymin:=y[1]- r[1];

d:=r[1]*2;

m:=0;

Randomize;

for i:=1 to n do {в цикле по числу историй}.

begin.

xp:=Random*d+xmin; {бросаем случайным образом}.

yp:=Random*d+ymin; {точки в выбранный квадрат}.

for j:=1 to 3 do {проверяем, попадает ли точка в каждый круг}.

if Sqr (xp-x[j])+Sqr (yp-y[j]) > rr[j] then goto NotInCircle;

Inc (m); {считаем количество точек, попавших сразу во все три круга}.

NotInCircle:

end;

WriteLn ('S = ', Sqr (d)*m/n:0:3); {результат, тем точнее, чем больше историй}.

end.

Задача 4. Два цилиндра одинакового радиуса R=1 пересекаются под прямым углом. Найти объем V их общей части.

Указание. В журнале «Квант» (№ 2, 1988 г.) приводится геометрический формализм решения задачи: V = 16R³/3.

Задача 5. Оценить чего больше: несократимых или сократимых дробей.

Фривольное условие задачи предполагает строгую формулировку: какова вероятность того, что наудачу взятая дробь несократима?

Ответ достаточно сложен и равен 6/р² = 0.6079… (Н.Я. Виленкин. В таинственном мире бесконечных рядов. Квант № 10, 1989).

Сформулируем условие иначе. Рассмотрим несократимые дроби вида a/b, где 1Ј a, b Ј.N. Количество их f (N). Нужно найти предел f (N) / N²для больших чисел N. Выберем случайные натуральные числа (не превосходящие фиксированного достаточно большого числа N) для числителя и знаменателя дроби. Повторяем «эксперимент» $n$ раз, подсчитывая количество (m) несократимых дробей, используя алгоритм Эвклида для нахождения наибольшего общего делителя числителя и знаменателя. Отношение m/$n$ дает оценку доли несократимых дробей.

Строго говоря, мы должны доказать, что характер соотношения не изменится при увеличении числа N. Мы же будем это предполагать во всех задачах, так как математические критерии, гарантирующие сходимость решения, известны в редких случаях. Сходимость можно проверять для различного числа испытаний (например, увеличивая в 10 раз). Важно, чтобы число историй по закону было «большим».

Упражнение 2.

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

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

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

Пример 3.

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

Указание. Пронумеруем клетки шахматной доски от 0 до 63. Случайное попадание зерна на какую-либо клетку моделируем случайным выбором клетки с помощью датчика случайных чисел RANDOM (64).

const.

Iterations = 10 000;

var.

Board: array[0.63] of LongInt;

Count: array[0.63] of LongInt;

i, j: LongInt;

begin.

Randomize;

FillChar (Count, SizeOf (Count), 0);

for i:=1 to Iterations do.

begin.

FillChar (Board, SizeOf (Board), 0);

for j:=0 to 63 do Inc (Board[Random (64)]);

for j:=0 to 63 do Inc (Count[Board[j]]);

end;

for j:=0 to 10 do.

WriteLn (j:2, ' ', Count[j]/Iterations:8:5);

end.

В результате моделирования оказывается, что вероятнее всего: 23 клетки останутся пустыми, на 24 клетках будет по одному зерну, на 12 клетках — по два, на 4 клетках — по три, на одной клетке — четыре.

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

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

На плоскости чертится прямоугольник, задающий границы лабиринта. Внутри прямоугольника выбирается точка (координаты которой задаются случайным образом), не лежащая на ранее построенных границах. От точки в случайном направлении (вправо, влево, вверх, вниз) рисуется линия границы до пересечения с какой-либо другой линией. Чтобы проходы в лабиринте были одинаковой ширины, координаты точки задаются с заранее выбранным шагом (например, на целочисленной сетке). Построение лабиринта прекращается по нажатию клавиши или когда выбраны все допустимые точки. Такой алгоритм построения не дает циклических путей в лабиринте и, следовательно, в нем всегда можно найти выход. На рис. 1 приведен вариант сгенерированного лабиринта.

Пример 4.

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

uses.

Graph, Crt;

const.

Step = 20; {шаг — размер клетки в пикселах}.

Width = 30; {ширина лабиринта в клетках}.

Height = 20; {высота лабиринта в клетках}.

dx: array[0.3] of Integer = (1, 0, -1, 0); {смещения по горизонтали}.

dy: array[0.3] of Integer = (0, -1, 0, 1); {и вертикали для четырех направлений}.

Рис. 1.

var.

Driver, Mode: Integer;

x, y, n: Integer;

Wall: Boolean;

begin.

Driver:=Vga; Mode:=VgaHi;

InitGraph (Driver, Mode, '');

{очертим границы лабиринта}.

Rectangle (0, 0, Width*Step, Height*Step);

Randomize;

repeat.

x:=Random (Width)*Step; {выбираем случайную точку}.

y:=Random (Height)*Step;

if GetPixel (x, y)0 then Continue; {если она лежит на стене, пробуем заново}.

MoveTo (x, y);

n:=Random (4); {выбираем случайное направление}.

repeat {рисуется линия от текущей точки}.

Inc (x, dx[n]*Step); {с заданным шагом и направлением}.

Inc (y, dy[n]*Step);

Wall:=(GetPixel (x, y)0); {до ближайшей стенки}.

LineTo (x, y);

until Wall;

until KeyPressed;

ReadKey;

CloseGraph;

end.