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

Интеграл как среднее

Вспомним: среднее значение функции \bar f = \frac{1}{b - a}\int_a^b f\,dx. Если X — случайная точка, равномерно распределённая на [a; b], то E f(X) = \bar f. По закону больших чисел среднее по многим случайным точкам приближается к математическому ожиданию:

Это метод Монте-Карло (назван в честь казино; придуман Уламом и фон Нейманом в Лос-Аламосе в 1940-х).

Метод попаданий

Площадь фигуры D внутри прямоугольника площади A: бросаем N случайных точек в прямоугольник, считаем, сколько попало в D:

Вычисление \pi: четверть круга в квадрате [0; 1]^2 имеет площадь \frac\pi4.

import numpy as np
import vis

rng = np.random.default_rng(0)
x, y = rng.random(2000), rng.random(2000)
inside = x**2 + y**2 <= 1
print("π ≈", 4 * inside.mean(), "по 2000 точкам")
vis.scatter(x, y, labels=inside.astype(int), title="Точки внутри четверти круга (оранжевые) и вне её", equal=True)

for N in [10**2, 10**3, 10**4, 10**5, 10**6]:
    p = (rng.random(N)**2 + rng.random(N)**2 <= 1).mean()
    print(f"N = {N:>8}: π ≈ {4 * p:.5f}, ошибка {abs(4 * p - np.pi):.5f}")

Ошибка \sim\frac{1}{\sqrt N}

По центральной предельной теореме ошибка оценки примерно нормальна со стандартным отклонением

где \sigma_f — стандартное отклонение f(X). Это значит:

• чтобы получить ещё один верный знак, нужно в 100 раз больше точек — в одномерном случае это очень медленно (Симпсон: в 10^{1/4} \approx 1{,}8 раза);
• ошибка не зависит от размерности! Для интеграла по 100-мерному кубу скорость та же.

Ошибку можно оценить по самой выборке: \sigma_f \approx s (выборочное отклонение). С вероятностью ≈95% истинное значение лежит в интервале \pm2\sigma_N.

import numpy as np
from math import gamma, pi

rng = np.random.default_rng(42)
# Объём единичного шара в размерности d: доля точек куба [−1, 1]^d внутри шара, умноженная на 2^d
for d in [2, 3, 5, 8, 10]:
    N = 400000
    pts = rng.uniform(-1, 1, size=(N, d))
    hit = (np.sum(pts**2, axis=1) <= 1)
    est = 2**d * hit.mean()
    se = 2**d * hit.std() / np.sqrt(N)
    exact = pi**(d / 2) / gamma(d / 2 + 1)
    print(f"d = {d:2}: V ≈ {est:.4f} ± {2 * se:.4f}   (точно {exact:.4f})")

С ростом размерности доля попаданий падает (шар занимает всё меньшую часть куба), и относительная ошибка растёт — сказывается \sigma_f. Это лечится умными способами выбора точек.

Как уменьшить дисперсию

• Выборка по значимости: точки бросают гуще там, где |f| велико, и делят на плотность: \int f\,dx = E\frac{f(Y)}{p(Y)}, Y \sim p. При p \approx \frac{|f|}{\int|f|} дисперсия почти нулевая.
• Стратификация: делим отрезок на части и в каждую бросаем поровну точек.
• Квазислучайные последовательности (Соболя, Холтона) заполняют пространство равномернее случайных; ошибка близка к \frac1N вместо \frac{1}{\sqrt N}.

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

Итог

• \int_a^b f\,dx \approx (b - a)\cdot\overline{f(X_i)} по случайным точкам; площади — долей попаданий.
• Ошибка \sim\frac{\sigma_f}{\sqrt N} и не зависит от размерности — Монте-Карло незаменим для многомерных интегралов.
• Дисперсию уменьшают выборкой по значимости, стратификацией и квазислучайными точками.