Метод Монте-Карло

Сколько квадратных километров воды в Капчагайском водохранилище? Обходить берег не нужно: возьмём снимок со спутника и будем бросать на него точки наугад. Доля точек, упавших в воду, умноженная на площадь кадра, и есть ответ — приблизительный, но с честно указанной погрешностью. Кадр в сцене — 3517,5 км², и с каждой новой точкой полоса вокруг оценки на графике сужается.

Посмотрите на график. После тысячи точек оценка — 1097 км² плюс-минус сотня. После миллиона — 1124 ± 3. Точек стало в тысячу раз больше, а погрешность уменьшилась в тридцать раз: √1000 ≈ 31,6. Это главный закон метода, и о нём вся статья. Вторая вещь на графике — пунктир «справочник» на 1847 км². Миллион точек к нему не приближается и не приблизится никогда. Почему — в разделе о границах метода; нажмите «Что считается водой» и посмотрите на восточный конец водохранилища.

Идея

Если искомое число — это вероятность или среднее значение чего-то случайного, его можно оценить опытом: разыграть случай n раз и взять долю удач или среднее результатов. Ошибка такой оценки убывает как σ/√n, где σ — разброс одного опыта. В формуле нет ни размерности задачи, ни сложности фигуры: только разброс и число опытов.

Площадь водохранилища — среднее значение функции «1, если точка в воде, иначе 0» по всему кадру, умноженное на площадь кадра. Ровно так же интеграл любой функции по области — среднее её значений в случайных точках, умноженное на объём области. Метод не требует ни формулы границы, ни её гладкости: только умения ответить, внутри точка или нет. Для снимка это цвет пикселя, для сложной детали — проверка по чертежу, для нейтрона в реакторе — розыгрыш его судьбы.

Почему 1/√n

Для точки «в воде или нет» σ² = p(1 − p), где p — истинная доля воды, у нас около 0,32. По центральной предельной теореме при большом n оценка распределена почти нормально, и в 95 % опытов она отклоняется от правды не больше чем на 1,96 · σ/√n. Отсюда привычный интервал p̂ ± 1,96 · √(p̂(1 − p̂)/n). Чтобы сузить его вдвое, нужно вчетверо больше точек; чтобы получить ещё один верный десятичный знак — в сто раз больше.

Медленно? Да. Но закон не зависит от того, в скольких измерениях лежит задача. В двух измерениях та же точность получается и сеткой — √n на √n точек. А в тридцати сетка с двумя точками вдоль каждой оси — уже миллиард точек, а одна точка вдоль оси — не сетка вовсе. Случайным точкам всё равно: σ/√n.

Честный ли интервал

«95 %» обещает вот что: если повторять весь опыт снова и снова, интервал накроет правду в 95 случаях из ста. Проверим обещание прямо: сто независимых опытов по n точек, у каждого свой интервал. Красные — те, что промахнулись мимо правды.

Для водохранилища при n = 100 привычная формула держит слово: промахов около пяти на сотню. Теперь выберите «пруд, 1 %» — водоём, который занимает сотую часть кадра, — и n = 30. Три опыта из четырёх не поймают ни одной точки в воде, и формула выдаст интервал нулевой ширины: «воды нет, и мы в этом уверены». Сколько это в точности, можно посчитать без всякой симуляции: число попаданий подчиняется биномиальному закону, и вероятность накрыть правду — сумма вероятностей тех k, чьи интервалы её накрывают.

   share      n    P(no hit)   textbook   Wilson
    lake    100        0.000      0.946    0.959
    pond     30        0.740      0.260    0.964
    pond    100        0.366      0.633    0.921

Обещано 95 %, а для пруда на 30 точках привычный интервал накрывает правду в 26 % случаев. В 2001 году Лоуренс Браун, Тони Цай и Анирбан ДасГупта показали, что покрытие этой формулы скачет неровно даже при больших n, а советы учебников, когда ей можно доверять, ненадёжны, и рекомендовали вместо неё интервал Уилсона.

Идея Эдвина Уилсона (1927) — не подставлять оценку p̂ в формулу разброса, а честно спросить: при каких значениях p наблюдённая доля p̂ правдоподобна, то есть отличается от p не больше чем на 1,96 разброса, посчитанного по самому p?

Возведите в квадрат — получится квадратное неравенство относительно p, и его решение — отрезок с концами

При p̂ = 0 он не схлопывается в точку, никогда не вылезает за [0, 1] и держит покрытие близко к обещанному: для того же пруда — 96,4 %. Переключите сцену на «Уилсон» и сравните.

Код

Программа читает маску воды — тот самый снимок, пересчитанный в равновеликую проекцию, где каждый пиксель покрывает одинаковые 97 × 97 м, — и бросает на него точки. Генератор случайных чисел свой и одинаковый в Python и Go, поэтому обе версии бросают одни и те же точки; первый дождь в сцене — тоже они.

frame 960x390 px = 3517.5 km2; water pixels 119639 = 1124.0 km2
        n   estimate, km2   ±95%      Wilson 95%         covers
       10          1055.2    999.1   [ 379.2, 2121.8]   yes
     1000          1097.4    101.0   [ 999.1, 1200.8]   yes
  1000000          1124.2      3.2   [1120.9, 1127.4]   yes

Десять точек дают «1055 ± 999» — честно, но бесполезно. Интервал Уилсона на тех же десяти точках, от 379 до 2122, несимметричен: он помнит, что площадь не бывает отрицательной.

Сколько это стоит

Время — O(n) испытаний, каждое стоит одну проверку «внутри или нет», память — O(1): хранить нужно только счётчик. Платят за точность: ошибка ∝ 1/√n, и каждый следующий знак стоит в сто раз больше. Миллион точек — секунды даже на Python; ещё два знака потребовали бы десяти миллиардов.

Этот закон можно обойти, если отказаться от полной независимости точек. Разрежем кадр на k × k клеток и бросим по одной случайной точке в каждую. Ошибка складывается теперь только из клеток, через которые проходит берег: клетки целиком в воде или целиком на суше дают точный ответ. Таких «береговых» клеток порядка k, то есть √n, и ошибка падает как n^−3/4 вместо n^−1/2 — этот результат для пикселей с краем есть у Дона Митчелла (1996), а сам приём в компьютерную графику принесли Роберт Кук, Томас Портер и Лорен Карпентер в 1984 году. Переключите первую сцену на «По клеткам» — на тех же n интервал заметно уже.

     n   plain: rms error, km2   stratified   ratio
   100                   139.7         85.4     1.6
  6400                    22.7          4.1     5.6

На 6400 точках клетки ошибаются в 5,6 раза меньше, и разрыв растёт с n. В многомерных задачах клетки строить трудно, и там используют их родственников — квазислучайные последовательности, которые раскладывают точки равномерно без сетки.

Где без него не обойтись

Площадь на плоскости можно посчитать и без случайности — сеткой или по формуле многоугольника. Метод Монте-Карло незаменим там, где измерений много. Вот задача с девяноста измерениями. Миллион тенге лежит в долларах, евро и юанях поровну. Сколько он будет стоить в тенге через месяц? Месяц — это 30 дней, и в каждый день меняются три курса: исход зависит от девяноста случайных чисел. Формулы нет: дневные изменения курсов не подчиняются ни одному аккуратному закону, а сетка с двумя точками на каждую из 90 осей — это 2⁹⁰ сценариев.

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

1899 daily moves, 20000 simulated months
mean after a month: 1,000,933 tenge
P(lose more than 5%) = 3.95% ± 0.27%
5% worst month: below 953,139 tenge, 95% interval [952,352, 954,076]
real 30-day windows losing more than 5%: 87 of 1870 = 4.65%

Вероятность потерять больше 5 % за месяц — 3,95 % ± 0,27 %. Самые плохие 5 % месяцев оставляют от миллиона меньше 953 тысяч; у этой границы тоже есть интервал: соседние по рангу исходы, от 952 до 954 тысяч. Так банки и оценивают рыночный риск портфеля — разыгрывают сценарии на истории.

А теперь последняя строка: в настоящей истории месяцев с такой потерей 4,65 % — больше, чем насчитала симуляция, и вне её интервала. Монте-Карло не ошибся в арифметике: он честно ответил на вопрос о модели, в которой дни независимы. Главная причина, по всей видимости, в том, что настоящие плохие месяцы — это серии плохих дней подряд, как весной 2022 года, когда доллар за месяц то дорожал на 20 %, то дешевел на 12 %, а перемешивание дней эти серии разбивает. Интервал описывает случайную ошибку и ничего не знает о том, верна ли модель.

Где метод бессилен

Систематическая ошибка. Вернёмся к водохранилищу. Миллион точек дал 1124 ± 3 км², а в справочнике — 1847 км². Метод здесь ни при чём: он сходится к тому, что считает водой классификатор, — к 1124,0 км² по всем пикселям снимка. Водой считается тёмный сине-зелёный пиксель, и мутный мелководный восточный конец, куда впадает река, на снимке 2016 года выглядит как суша. А к какому уровню воды относятся справочные 1847 км², справочник не говорит; известно, что водохранилище так и не наполнили до проектного объёма. Сколько ни добавляй точек, такую ошибку не уменьшить: она в определении, а не в случайности.

Редкие события. Чтобы оценить вероятность 10⁻⁶ с относительной ошибкой 10 %, нужно, чтобы событие случилось хотя бы сотню раз, — сто миллионов опытов. Относительная ошибка оценки маленькой вероятности p — примерно 1/√(np): для p = 10⁻¹² метод бессилен. Такие вероятности считают точно, как во второй задаче, или меняют распределение так, чтобы редкое стало частым, и пересчитывают веса — это выборка по значимости.

Плохая случайность. Всё сказанное держится на независимости испытаний. Генератор с короткой или узорчатой последовательностью ломает оценку незаметно; о том, какие генераторы годятся, — статья о генераторах случайных чисел.

Где это работает

• Кино и игры. Трассировка путей считает освещённость каждого пикселя как среднее по случайным путям света, и шум на недоделанном кадре — это ошибка Монте-Карло: чтобы уменьшить его вдвое, лучей нужно вчетверо больше. Так работает, например, Cycles — движок рендеринга Blender. Разбиение пикселя на клетки с одним случайным лучом в каждой — та самая стратификация.
• Финансы. Риск портфеля, цена сложных опционов, сценарии платежей — всё, где исход зависит от многих случайных величин сразу.
• Физика. Перенос нейтронов и фотонов, молекулы, рассеяние в детекторах — с этой задачи метод и начался.
• Геометрия и инженерия. Объёмы и площади тел, которые заданы только проверкой «внутри или снаружи», допуски сборки из многих деталей.

История

Идея пришла Станиславу Уламу в 1946 году за пасьянсом: вероятность, что он сойдётся, проще оценить, разложив его сотню раз, чем посчитать комбинаторикой. Вместе с Джоном фон Нейманом он применил её к переносу нейтронов, а название по казино предложил их коллега Николас Метрополис. Первая статья с названием «Метод Монте-Карло» — Метрополиса и Улама, 1949 год. Подробно эта история рассказана в главе курса о рандомизированных алгоритмах, а там же есть и самый известный опыт метода — число π дротиками. Интервал для доли, которым мы пользовались, Эдвин Уилсон опубликовал в 1927 году — за двадцать лет до того, как у метода появилось имя.

Задачи

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