Метод Монте-Карло: интеграл через случайность

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

Метод Монте-Карло оценивает интеграл как среднее значение функции в случайных точках, умноженное на объём области: ∫f dV ≈ V·(среднее f по случайным точкам).

Идея: средним по случайным точкам

Интеграл функции по области — это, по сути, её среднее значение, умноженное на объём области. А среднее можно оценить, набросав случайные точки и усреднив значения функции в них (закон больших чисел). Чем больше бросков, тем точнее оценка. Звучит грубо — и в одномерном случае Монте-Карло действительно проигрывает трапециям. Но у него есть суперсила, ради которой его и придумали (для атомной бомбы в Лос-Аламосе): он почти не зависит от размерности.

Классическая демонстрация — оценка числа π: бросаем точки в единичный квадрат и считаем долю попавших в четверть круга. Эта доля — отношение площадей, то есть π/4.

import random

def оценка_пи(n):
    random.seed(42)
    внутри = 0
    for _ in range(n):
        x = random.random()
        y = random.random()
        if x*x + y*y <= 1.0:    # попала в четверть круга радиуса 1
            внутри += 1
    return 4.0 * внутри / n

import math
print("оценка π методом Монте-Карло:")
for n in [100, 10000, 1000000]:
    оценка = оценка_пи(n)
    print(f"  n={n:>7}: π ≈ {оценка:.5f}  (ошибка {abs(оценка - math.pi):.5f})")

Вывод:

оценка π методом Монте-Карло:
  n=    100: π ≈ 3.12000  (ошибка 0.02159)
  n=  10000: π ≈ 3.12600  (ошибка 0.01559)
  n=1000000: π ≈ 3.14059  (ошибка 0.00100)

Оценка приближается к π, но сходится медленно — в среднем как 1/√n. От n=100 до n=10⁶ точность улучшилась лишь примерно в 20 раз, хотя точек стало в 10000 раз больше: закон 1/√n жесток. Монте-Карло к тому же вероятностен — отдельный прогон может «повезти» или «не повезти», поэтому оценки пляшут вокруг истины, а не убывают строго монотонно.

Скорость сходимости: 1/√n

Ошибка Монте-Карло убывает как 1/√n — независимо от размерности. Это медленно: для одной дополнительной цифры (ошибка ÷10) нужно n × 100. В 1D трапеции с их h² бьют Монте-Карло наголову. Но вот ключевой момент: у сеточных методов число узлов растёт как m^d (где d — размерность), и в 10-мерном интеграле даже 10 узлов на ось дают 10^10 точек — «проклятие размерности». А ошибка Монте-Карло остаётся 1/√n при любом d. Поэтому для интегралов размерности выше ~4 Монте-Карло часто единственный реальный путь.

Как работает под капотом

Точность Монте-Карло определяется дисперсией функции: ошибка ≈ σ/√n, где σ — стандартное отклонение значений f. Отсюда главные приёмы ускорения — уменьшение дисперсии: importance sampling (брать точки чаще там, где f велика), стратификация (делить область на ячейки), квази-Монте-Карло (детерминированные «равномерно рассыпанные» последовательности Соболя вместо чистого random — они дают сходимость почти 1/n). В машинном обучении и финансах (оценка опционов, байесовский вывод) интегралы бывают в сотнях измерений — там Монте-Карло и его варианты безальтернативны.

Частые ошибки

• Применять Монте-Карло к гладкому 1D-интегралу. Там трапеции/Симпсон точнее на порядки; Монте-Карло оправдан в высоких размерностях.
• Ждать точности без огромного n. Сходимость 1/√n медленна; для высокой точности нужны миллионы–миллиарды точек или снижение дисперсии.
• Не фиксировать seed при отладке. Случайность мешает воспроизводимости; для тестов задавайте random.seed.

Итоги

• Монте-Карло оценивает интеграл средним значением функции по случайным точкам, умноженным на объём области.
• Сходимость 1/√n — медленная, но не зависит от размерности.
• В 1D проигрывает сеточным методам; в высоких размерностях (проклятие размерности) — часто единственный путь.
• Ускоряют снижением дисперсии: importance sampling, стратификация, квази-Монте-Карло (Соболь).