Оценка интегралов методом Монте-Карло

Применяем случайность к вычислению интегралов — в том числе там, где обычные методы пасуют.

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

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

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

Рассмотрим интеграл на отрезке $[a;b]$. Если брать точки $x$ равномерно на этом отрезке, то среднее значение функции $f(x)$ умноженное на длину отрезка даёт интеграл:

$$\int_a^b f(x)\,dx=(b-a)\,\mathbb{E}[f(X)]\approx \frac{b-a}{n}\sum_{i=1}^{n} f(x_i).$$

Идея в том, что $\mathbb{E}[f(X)]$ для равномерной $X$ — это среднее высоты графика, а площадь прямоугольника «средняя высота на ширину» и есть интеграл. Проверим на интеграле, который легко взять точно:

$$\int_0^1 x^2\,dx=\frac{1}{3}\approx 0{,}3333.$$

import random
random.seed(51)

def f(x):
    return x * x

n = 2000000
total = sum(f(random.random()) for _ in range(n))
estimate = (1 - 0) * total / n
print("Монте-Карло:", round(estimate, 4))
print("Точно 1/3: ", round(1/3, 4))

Вывод:

Монте-Карло: 0.3334
Точно 1/3:  0.3333

Оценка сошлась к $\frac{1}{3}$. Мы вычислили интеграл, ни разу не взяв первообразную, — только усреднив значения функции.

Когда метод незаменим

Возьмём интеграл, у которого нет первообразной в элементарных функциях:

$$\int_0^1 e^{-x^2}\,dx\approx 0{,}7468.$$

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

import random, math
random.seed(52)

n = 2000000
total = sum(math.exp(-random.random()**2) for _ in range(n))
print("Монте-Карло:", round(total / n, 4))
print("Эталон:    ", 0.7468)

Вывод:

Монте-Карло: 0.7469
Эталон:     0.7468

Оценка совпала с табличным значением до четвёртого знака. В этом сила метода: ему всё равно, есть ли у функции первообразная.

Сила в многомерности

На прямой Монте-Карло уступает обычным численным методам (метод прямоугольников, Симпсона), которые точнее при том же числе точек. Но с ростом размерности классические сетки взрываются: для $d$ измерений нужно $k^d$ узлов. Ошибка же Монте-Карло убывает как $1/\sqrt{n}$ независимо от размерности — поэтому в задачах с десятками переменных он остаётся практически единственным выходом.

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

Код генерирует случайные точки на $[0;1]$ и усредняет значения функции в них — это прямая численная реализация $\mathbb{E}[f(X)]$. Умножение на длину отрезка превращает среднюю высоту в площадь. По закону больших чисел выборочное среднее $\frac{1}{n}\sum f(x_i)$ сходится к истинному ожиданию, а значит, к интегралу. Метод снова оказывается переодетым законом больших чисел: интеграл — это ожидание, а ожидание мы оцениваем средним.

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

Первая ошибка — забыть множитель $(b-a)$: на отрезке $[0;1]$ он равен 1 и «прячется», но на других отрезках без него ответ будет в разы меньше. Вторая — применять формулу с равномерными точками к неравномерной выборке: тогда нужно делить на плотность (метод важности), иначе оценка смещается. Третья — ждать высокой точности на одномерных гладких функциях: тут обычные методы лучше, Монте-Карло раскрывается в больших размерностях.

Итог

Интеграл равен длине области, умноженной на среднее значение функции.
Монте-Карло берёт интегралы без первообразной, например $\int_0^1 e^{-x^2}dx$.
Ошибка убывает как $1/\sqrt{n}$ независимо от размерности.
Метод незаменим в задачах высокой размерности, где сетки взрываются.