Интегрирование методом Монте-Карло

Определённый интеграл — это, по сути, среднее значение функции, растянутое на длину отрезка. А среднее прекрасно оценивается случайными выстрелами.

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

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

Интеграл как среднее значение

Вспомним смысл определённого интеграла. Величина ∫f(x)dx на отрезке от a до b — это площадь под графиком функции на этом отрезке. С другой стороны, есть понятие среднего значения функции на отрезке: если бы график был «выровнен» в горизонтальную прямую, на какой высоте она бы прошла, сохранив ту же площадь? Связь между ними прямая:

интеграл = (среднее значение функции на отрезке) × (длина отрезка)

Эта формула — ключ ко всему. Длину отрезка (b − a) мы знаем точно. А среднее значение функции легко оценить методом Монте-Карло: выберем много случайных точек на отрезке, посчитаем в них значения функции и возьмём арифметическое среднее. По закону больших чисел это среднее сходится к истинному среднему значению функции.

Считаем интеграл x² на отрезке [0, 1]

Проверим метод на примере, где ответ известен заранее. Интеграл функции x² от 0 до 1 равен 1/3 ≈ 0.33333. Длина отрезка здесь равна 1, поэтому интеграл совпадает со средним значением функции — нам достаточно усреднить x*x по случайным точкам.

import random
random.seed(1)
N = 50000
s = 0.0
for _ in range(N):
    x = random.random()
    s += x * x
print(f"Оценка интеграла x^2 на [0,1]: {s/N:.5f}")
print("Точное значение: 0.33333")

Вывод:

Оценка интеграла x^2 на [0,1]: 0.33396
Точное значение: 0.33333

Полсотни тысяч случайных точек дали третий знак почти верно. В переменной s мы накапливаем сумму значений x*x, а s/N — это и есть среднее. Поскольку длина отрезка равна 1, отдельно домножать ни на что не нужно. Если бы отрезок был, скажем, от 0 до 2, мы бы генерировали x = random.random() * 2 и в конце умножали среднее на длину 2.

Общая схема для любого отрезка

Чтобы метод работал на произвольном отрезке от a до b, достаточно двух правок: растягивать случайную точку на нужный диапазон и умножать среднее на длину отрезка.

import random
random.seed(7)

def integrate(f, a, b, N=50000):
    total = 0.0
    for _ in range(N):
        x = a + (b - a) * random.random()
        total += f(x)
    return (b - a) * total / N

print(f"x^2 на [0,2]: {integrate(lambda x: x*x, 0, 2):.4f}")
print("Точное значение: 2.6667")

Вывод:

x^2 на [0,2]: 2.6579
Точное значение: 2.6667

Где метод по-настоящему незаменим: многомерные интегралы

Для одномерного интеграла, честно говоря, есть способы и поточнее — метод прямоугольников, трапеций, Симпсона. Они сходятся к ответу быстрее, чем Монте-Карло. Так зачем тогда вообще случайность?

Ответ — в размерности. Представьте интеграл не по отрезку, а по объёму в пространстве многих измерений: десять, сто, тысяча переменных сразу. Такие интегралы реально возникают в физике, финансах и машинном обучении. Классические методы здесь беспомощны: чтобы покрыть сеткой пространство из 100 измерений, понадобилось бы астрономическое число узлов — это называют «проклятием размерности». А Монте-Карло почти не замечает числа измерений: формула «среднее значение функции × объём области» работает одинаково и для отрезка, и для стомерного куба. Достаточно уметь генерировать случайную точку в этом пространстве и считать в ней функцию. Именно поэтому в серьёзных расчётах Монте-Карло часто оказывается единственным реально работающим инструментом.

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

Разберём цепочку по шагам. Мы хотим оценить I = ∫f(x)dx на отрезке от a до b. Определение среднего значения функции даёт среднее = I / (b − a), откуда I = (b − a) × среднее. Среднее значение функции по отрезку — это математическое ожидание величины f(x), где x равномерно распределена на отрезке. А математическое ожидание Монте-Карло оценивает арифметическим средним по выборке.

x1, x2, ..., x_N  — равномерно случайны на [a, b]
        |
        v
среднее ≈ (f(x1) + f(x2) + ... + f(x_N)) / N
        |
        v
интеграл ≈ (b − a) × среднее

В коде это ровно три действия: накопить сумму total значений функции, поделить на N (получить среднее), умножить на длину (b − a). Когда длина равна 1, последний шаг невидим — поэтому в первом примере его и не было.

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

Забыть умножить среднее на длину отрезка (b − a). На отрезке [0, 1] ошибка незаметна (длина равна 1), но на любом другом отрезке ответ окажется в несколько раз меньше или больше.
Генерировать точки не на том диапазоне: использовать random.random() (это всегда отрезок от 0 до 1) там, где нужен другой отрезок, не растянув точку формулой a + (b − a) * random.random().
Применять Монте-Карло для простого одномерного интеграла и ждать высокой точности — там методы трапеций и Симпсона выгоднее. Сила Монте-Карло проявляется в многомерье.
Делить на неверное число: усреднять надо именно по числу испытаний N, а не по чему-то ещё.
Накапливать сумму в целочисленной переменной. Начинайте с 0.0, иначе при делении легко потерять дробную часть при неосторожном коде.

Итоги

Определённый интеграл равен среднему значению функции, умноженному на длину области интегрирования.
Монте-Карло оценивает это среднее как арифметическое среднее значений функции в случайных точках.
На отрезке [0, 1] длина равна 1, поэтому интеграл совпадает со средним и домножение не требуется.
Для произвольного отрезка точку растягивают формулой a + (b − a) * random.random(), а среднее умножают на (b − a).
Главная сила метода — многомерные интегралы, где классические сеточные методы разбиваются о «проклятие размерности».