Каждый знает, как рассчитать число π методом Монте-Карло. Об этом упоминают ещё в школе на уроке информатики. Но за скобками обычно остаётся один небольшой вопрос, о котором будет сказано чуть дальше.
Позволю себе напомнить суть метода Монте-Карло применительно к вычислению числа π "школьным" способом. Берётся единичная окружность (окружность с центром в начале координат и радиусом, равным единице) и вписывается в квадрат. Сторона у этого квадрата получается равной двум.
Затем из кармана вытаскивается датчик случайных чисел, и с помощью этого датчика генерируется N точек, равномерно разбросанных внутри квадрата. Наконец, подсчитывается количество M точек, попавших внутрь окружности, и отношение M : N дает приблизительную оценку отношения площади круга (равной π) к площади вмещающего его квадрата (равной 4).
Маленькое замечание относительно получения этих случайных точек. Как правило, компьютер предоставляет (или эмулирует) так называемый датчик случайных чисел, который выдаёт последовательность чисел - реализаций случайной величины, равномерно распределенной на интервале [0, 1). С помощью такого датчика можно получить координаты очередной точки по простым формулам:
x = 2 ∙ <одно_значение_датчика_сл_ч> − 1
y = 2 ∙ <ещё_одно_значение_датчика_сл_ч> − 1
Хорошо, у нас есть алгоритм, который постепенно сходится к значению π. Теперь настало время задаться тем самым вопросом, о котором говорилось вначале: сколько точек потребуется накидать в квадрат, чтобы получилось значение π с интересующей нас точностью?
Тут надо сделать небольшое лирическое отступление. Алгоритмы метода Монте-Карло, в отличие от прочих, имеют вероятностный характер, то есть, если нам очень сильно не повезёт, то все наши точки окажутся в круге, и мы получим значение π, равное четырём. Или наоборот, в круг не попадёт ни одна из N разыгранных точек, и на выходе окажется π, равное нулю. Поэтому задача оценки зависимости числа N необходимых испытаний от точности может быть поставлена лишь следующим образом:
Каким должно быть число испытаний N, чтобы рассчитанное значение лежало в пределах заданной окрестности (π−ε, π+ε) истинного значения π с заданной вероятностью 1−δ ? Например, сколько точек нам понадобится, чтобы получилось число π, верное хотя бы до второго знака (т.е., ε = 0.05), с вероятностью 90% (т.е., δ = 0.1) ?
Ответ на этот вопрос потребует экскурса в теорию вероятностей. Рассмотрим случайную величину ξ, принимающую значение 1, если случайная точка попала в круг, и 0 в противном случае. Эта случайная величина имеет так называемое распределение Бернулли, задаваемое вероятностями:
P { ξ = 1 } = p = π / 4
P { ξ = 0 } = q = 1 − p = 1 − π / 4
У этой случайной величины есть пара важных характеристик. Они называются математическое ожидание и дисперсия (семантикой этих слов на данном этапе можно пренебречь). Математическое ожидание (обозначается Eξ) и дисперсия (обозначается Dξ) выражаются следующим образом:
Eξ = 0 ∙ P { ξ = 0 } + 1 ∙ P { ξ = 1 } = 0 ∙ q + 1 ∙ p = p
Dξ = (0 − Eξ)² ∙ P { ξ = 0 } + (1 − Eξ)² ∙ P { ξ = 1 } = p² ∙ q + q² ∙ p = p ∙ q
Обозначим ζ cумму из N таких случайных величин. Эта новая случайная величина имеет так назваемое биномиальное распределение, и её математическое ожидание и дисперсия выражаются так:
Eζ = N ∙ p
Dζ = N ∙ p ∙ q
Очевидно, наше бросание N случайных точек в квадрат с подсчётом количества точек, попавших в круг, является реализацией случайной величины ζ' = ζ / N. Её математическое ожидание и дисперсия, в свою очередь, выглядят так:
Eζ' = p
Dζ' = p ∙ q / N
А дальше воспользуемся знаменитым неравенством Чебышёва, которое гласит:
P { | ζ' − Eζ' | > ε } < Dζ' / ε²
Или, говоря человеческим языком, вероятность того, что реализация случайной величины отклонится от своего математического ожидания больше, чем на ε, оказывается меньше, чем дисперсия этой величины, делённая на квадрат этого ε.
"Перевернув" это неравенство получим:
P { | ζ' − Eζ' | < ε } > 1 − Dζ' / ε²
Теперь приравняем правую часть к заданному δ и получим уравнение для определения того, каким должно быть искомое число точек N:
Dζ' / ε² ≤ δ
p ∙ q / N ≤ δ ∙ ε²
N ≥ p ∙ q / (δ ∙ ε²)
В нашем случае, учитывая, что реализация случайной величины ζ' дает оценку не π, а π / 4, получаем:
N ≥ π / 4 ∙ ( 1 − π / 4 ) / ( δ ∙ (ε / 4)² ) = π ∙ ( 4 − π ) / ( δ ∙ ε² )
Так как величина π нам заранее неизвестна, можно вместо неё использовать "грубую" оценку 3 < π и получить искомую оценку для необходимого числа точек N:
N ≥ 3 / ( δ ∙ ε² )
Или, в числах: N ≥ 3 / ( 0.1 ∙ 0.05² ) = 12'000 точек.
Это, конечно, немало, но основная печаль начинается дальше. Если мы захотим получить ещё один верный знак, то число ε надо будет уменьшить в 10 раз, и, соответственно, число N возрастёт в сто(!) раз.
В этом заключена и слабость метода Монте-Карло, и его сила. Дело в том, что подобное поведение метода - зависимость погрешности от числа испытаний как ε ~ 1/√N - практически не зависит ни от сложности вычисляемого отношения площадей, ни от размерности задачи. То есть, метод Монте-Карло сходится медленно - но верно.
Ну и, конечно, нельзя не сказать, что вышеупомянутая задача - это даже не вершина айсберга, а... так, маленький камешек на берегу материка под названием "статистическое моделирование".