Оценка числа π методом Монте-Карло

Как с помощью случайных точек и интеграла оценить число π прямо в браузере.

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

Иногда честно посчитать интеграл сложно, особенно в многомерье. Тогда выручает неожиданно простой приём — бросать случайные точки и считать, какая доля попала «внутрь». Это метод Монте-Карло, и он позволяет элегантно оценить число $\pi$. Идея в том, что площадь четверти круга радиуса $1$ равна $\frac{\pi}{4}$, и эту площадь можно записать интегралом.

Связь π с интегралом

Площадь четверти единичного круга — это интеграл функции верхней дуги:

$$\int_0^1 \sqrt{1-x^2}\,dx = \frac{\pi}{4}$$

Значит, $\pi = 4\int_0^1 \sqrt{1-x^2}\,dx$. Можно посчитать этот интеграл суммой Римана, а можно — методом Монте-Карло, который особенно нагляден.

Монте-Карло: бросаем точки

Кидаем случайные точки в единичный квадрат $[0,1]\times[0,1]$. Точка попадает в четверть круга, если $x^2+y^2 \le 1$. Доля попавших точек примерно равна отношению площадей — то есть $\frac{\pi}{4}$. Умножаем долю на $4$ и получаем оценку $\pi$.

import random
import math

random.seed(42)
for n in [1000, 10000, 100000, 1000000]:
    inside = 0
    for _ in range(n):
        x = random.random()
        y = random.random()
        if x*x + y*y <= 1.0:
            inside += 1
    pi_est = 4 * inside / n
    print(f"n={n:>7}  оценка pi = {pi_est:.5f}  ошибка={abs(pi_est-math.pi):.5f}")

print("настоящее pi =", round(math.pi, 5))

Вывод:

n=   1000  оценка pi = 3.12800  ошибка=0.01359
n=  10000  оценка pi = 3.12880  ошибка=0.01279
n= 100000  оценка pi = 3.13988  ошибка=0.00171
n=1000000  оценка pi = 3.14056  ошибка=0.00103
настоящее pi = 3.14159

С ростом числа бросков оценка приближается к $\pi$. Случайность, усреднённая по миллиону точек, выдаёт фундаментальную математическую константу. Зерно генератора фиксировано (`seed(42)`), поэтому вывод воспроизводим.

Скорость сходимости

У Монте-Карло своя цена: ошибка убывает медленно — как $\frac{1}{\sqrt{n}}$. Чтобы добавить один верный знак, надо в $100$ раз больше точек. Для одномерных интегралов это хуже трапеций и Симпсона. Но в высоких размерностях Монте-Карло незаменим: там сетка из прямоугольников взрывается экспоненциально, а случайные точки работают всегда.

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

Каждая точка — это честный «эксперимент»: попала в круг или нет. Доля успехов — это эмпирическая вероятность попадания, а вероятность по определению равна отношению благоприятной площади ко всей. Поскольку площадь квадрата равна $1$, доля попаданий прямо даёт площадь четверти круга. Закон больших чисел гарантирует, что с ростом числа испытаний эмпирическая доля сходится к истинной вероятности.

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

Первая — забыть умножить долю на $4$: сама доля даёт лишь $\frac{\pi}{4}$. Вторая — ждать высокой точности при малом $n$: из-за закона $\frac{1}{\sqrt{n}}$ тысяча точек даёт лишь два-три знака. Третья — не зафиксировать зерно генератора и удивляться, что вывод каждый раз разный (для воспроизводимости ставьте `random.seed`). Четвёртая — применять Монте-Карло там, где детерминированные методы (Симпсон) точнее и быстрее.

Итог

• Монте-Карло оценивает площади и интегралы случайными точками.
• $\pi = 4\cdot$ доля точек, попавших в четверть круга.
• Ошибка убывает как $\frac{1}{\sqrt{n}}$ — медленно, но надёжно.
• Метод незаменим в высоких размерностях.