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

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

Метод Монте-Карло оценивает интеграл как среднее значение функции в случайных точках, умноженное на объём области: $\int f \, dV \approx V \cdot (\text{среднее } f \text{ по случайным точкам})$.

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

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

Классическая демонстрация — оценка числа $\pi$: бросаем точки в единичный квадрат и считаем долю попавших в четверть круга. Эта доля — отношение площадей, то есть $\pi/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)

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

Скорость сходимости: $1/\sqrt{n}$

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

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

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

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

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

Итоги

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