Материал: Имитационное и статистическое моделирование — копия

Формула (1.41), опирающаяся на значение дисперсии D(x) усредняемой с. в. x, является основной формулой планирования статистических экспериментов.

При вычислении эмпирической вероятности усредняемая с. в. x является

известна: p = S = π/4 = 0,785398, откуда D(x) = p(1–p) = 0,168548. Поскольку n = 10 000, то из (1.41) получаем σn ≈ 0,004, 3σn ≈ 0,012. Нажимая многократно клавишу F9 для пересчёта листа, можно убедиться, что абсолютная величина погрешности в ячейке F1 действительно практически никогда не превышает

0,012.

Следует заметить, что на практике, когда применяется метод Монте-Карло, точное значение вероятности p, конечно, не бывает известно. Поэтому вероятность p при расчёте D(x) определяется приближённо, уже в ходе применения метода. В рассмотренном только что случае эта приближённая вероятность находится в ячейке D1. Её использование вместо точной вероятности не изменяет существенно выводов о погрешности расчёта.

Пример 1 2 . Сравним погрешность рассмотренного варианта метода Монте-Карло с погрешностью пошаговых методов интегрирования. Поскольку погрешность интегрирования зависит от числа точек в области интегрирования, по которым проводится вычисление, сравнение будем производить при одинаковом числе точек n = 10 000. Расчёт интеграла как вероятности будем сравнивать с детерминированным расчётом методом сеток [26].

Вначале сравним расчёт площади S (рис. 1.5), уже проведённый методом Монте-Карло, с расчётом этой же площади методом сеток. На чистом листе Ms Excel введём координаты 10 000 точек, заполняющих единичный квадрат

спостоянным шагом по вертикали и по горизонтали, равным 0,01. Квадрат должен быть заполнен 100 рядами точек по 100 точек в каждом ряду. Наиболее быстро это можно сделать следующим образом.

Вдиапазон A1:A10000 занесите арифметическую прогрессию от 0 до 9999

сшагом 1 (для этого запишите число 0 в ячейку A1 и используйте меню Правка/Заполнить/Прогрессия…). Далее используем, с небольшой модификацией, стандартный программистский приём преобразования адресов элементов линейного массива в адреса элементов двумерного массива. В ячейку B1 запишите формулу =ЦЕЛОЕ(A1/100)/100, а в ячейку C1 – формулу =ОСТАТ(A1;100)/100. Скопируйте B1 и C1 вниз на все 10 000 строк. Таким образом, теперь в диапазоне B1:C10000 построчно записаны координаты всех 10 000 точек, регулярной сеткой покрывающих единичный квадрат. Эту сетку можно увидеть, построив по диапазону B1:C10000 диаграмму типа Точечная. Чтобы маркеры на ней не сливались, приходится уменьшать размер маркеров до 2 пт, растягивать диаграмму, и/или увеличивать масштаб отображения листа.

Для вычисления интеграла в ячейку D1 введите формулу, определяющую, находится ли точка регулярной сетки в круге: =ЕСЛИ(B1^2+C1^2<1;1;0),

26

и скопируйте её вниз методом двойного щелчка. В ячейку F1 введите формулу =СУММ(D1:D10000)/10000. Здесь получен результат расчёта интеграла методом сеток. Он составляет 0,7949 и отличается от точного значения на 0,0095. Погрешность метода Монте-Карло, применённого нами выше, составила 0,012. Таким образом, оба сравниваемых метода характеризуются здесь сопоставимыми погрешностями.

Теперь сравним возможности расчёта этими методами кратных интегралов, решая задачу, подобную только что рассмотренной. В четырёхмерном пространстве положительный гипероктант содержит 1/16 часть четырёхмерного шара единичного радиуса, с центром в начале координат. Объём V этой части гипершара целиком заключён внутри единичного гиперкуба.

Расчёт объёма V шестнадцатой части гипершара методом Монте-Карло почти не отличается от выполненного выше расчёта площади S четверти круга. Отличия состоят лишь в том, что четырьмя координатами случайных точек =СЛЧИС() заполняются четыре столбца A1:D10000, а формула индикатора попадания точки в гипершар, которую можно записать в ячейку E1, имеет вид =ЕСЛИ(A1^2+B1^2+C1^2+D1^2<1;1;0). После её копирования в столбец E1:E10000 и расчёта объёма V как среднего значения в этом столбце получаем оценку объёма, составляющую примерно 0,31. Оценка погрешности по правилу трёх сигм составляет 0,014.

Теперь проведём расчёт объёма V методом сеток. Десять тысяч точек в регулярной сетке, заполняющей единичный гиперкуб, должны следовать друг за другом вдоль каждого из четырёх координатных направлений с шагом 0,1. На чистом листе Excel в диапазон A1:A10000 запишите числа от 0 до 9999 – это номера точек. Их четыре координаты найдите в столбцах B1:E10000 модифицированным стандартным приёмом программистов, использованным выше для двухмерного случая, а именно введите в ячейки B1, C1, D1 и E1 следующие четыре формулы:

=ЦЕЛОЕ(A1/1000)/10 =ЦЕЛОЕ(ОСТАТ(A1;1000)/100)/10 =ЦЕЛОЕ(ОСТАТ(A1;100)/10)/10 =ОСТАТ(A1;10)/10

В ячейку F1 введите формулу, определяющую, находится ли точка регулярной сетки внутри гипершара: =ЕСЛИ(B1^2+C1^2+D1^2+E1^2<1;1;0). Скопируйте все пять формул ячеек B1:F1 вниз методом двойного щелчка. В поле G1 результат расчёта определите по формуле =СУММ(F1:F10000)/10000. Он состав-

ляет 0,4213.

Итак, метод Монте-Карло дал оценку V ≈ 0,31 ± 0,014, пошаговый метод сеток при тех же затратах компьютерного времени (при том же количестве вычислений n = 10 000) показывает результат V ≈ 0,4213. Поскольку здесь вычисляется несложный интеграл, можно найти и точное его значение. Определяя посредством точного интегрирования объём четырёхмерного гипершара радиуса R, составляющий π2R4/2, находим, что при R = 1 его шестнадцатая часть V = π2R4/32 = 0,308425. Таким образом, при вычислении кратного интеграла метод Монте-Карло оказался более точным. Погрешность метода сеток (0,4213–0,308425) ≈ 0,113 оказалась приблизительно на порядок больше погрешности метода Монте-Карло.

27

Чтобы завершить характеристику достоинств метода Монте-Карло, предположим ещё, что погрешность выполненного расчёта необходимо уменьшить в k раз. При использовании метода Монте-Карло погрешность убывает примерно пропорционально квадратному корню из числа опытов n, поэтому n придётся увеличить в k2 раз. При использовании метода сеток погрешность убывает примерно пропорционально шагу между точками сетки, поэтому число точек придётся увеличить в k4 раз. Например, чтобы в рассмотренном примере расчёта пошаговый метод сеток достиг точности метода Монте-Карло, т. е. чтобы снизить погрешность метода сеток на порядок, нужно число точек в регулярной сетке увеличить приблизительно до n = 10 000×104 = 100 000 000.

с конечным или бесконечным интервалом интегрирования (a, b). Чтобы применить метод Монте-Карло, перепишем интеграл (1.42) в форме (1.33) м. о. функции и рассчитаем это м. о., используя оценку (1.37).

Вначале введём под знак интеграла какую-либо (произвольную) п. р. в. f(t):

Очевидно, выбранная функция f(t) не должна быть равна нулю в точках интервала интегрирования (a, b). Теперь для приведения интеграла к форме (1.33) остаётся, не меняя значения интеграла, перейти к бесконечным пределам интегрирования.

Для этого определим вспомогательную функцию H(t) через заданную подынтегральную функцию h(t) так, чтобы H(t) совпадала с h(t) на интервале (a, b), но была равна нулю вне этого интервала:

Тогда преобразование (1.43) можно продолжить так:

Представление (1.45) интеграла (1.42) является м. о. функции [H(x)/f(x)] от с. в. x ~ f(t) (сравните (1.45) с формулой (1.33)). Таким образом, искомый интеграл

28

т. е. может быть рассчитан методом Монте-Карло с использованием оценки (1.37). Для такого расчёта нужно только сгенерировать выборку с. в. x, имеющей п. р. в. f(t), вычислить и усреднить функцию [H(x)/f(x)] выборочных значений x.

Заметим, что если выбирается такая f(t), которая равна нулю при всех t (a,b) , то в представлении (1.46) вместо H(x) можно использовать непосредственно функцию h(x), поскольку в этом случае в генерируемой выборке x всегда попадает в (a, b).

4.4. Примеры расчета интеграла как математического ожидания

Пример 1 3 . Представим интеграл (1.39) в форме м. о., выбирая п. р. в. f(t) равномерного распределения R[0, 1], т. е. полагая f(t) = 1 при 0 ≤ t ≤ 1. Поскольку интервал распределения здесь совпадает с интервалом интегрирования, то нет необходимости использовать вспомогательную функцию H(x), и (1.46) записывается в виде

На чистом листе Ms Excel в диапазоне A1:A10000 сгенерируйте выборку

значений БСВ x. В ячейке B1 реализуйте усредняемую в (1.47) с. в. 1 − x2 , используя для этого формулу =КОРЕНЬ(1-A1^2). Затем заполните выборкой усредняемых значений с. в. диапазон B1:B10000, копируя ячейку B1 методом двойного щелчка.

Вячейке D1 вычислите оценку =СУММ(B1:B10000)/10000 интеграла S,

вячейке D2 – оценку =СТАНДОТКЛОН(B1:B10000) среднеквадратичного отклонения усредняемой с. в., необходимую для контроля точности. Фактическую погрешность расчёта =D1-ПИ()/4 введите в ячейку F1, а формулу предельного размаха погрешности =3*D2/100 введите в ячейку F2. Полезно ограничить представление числа в ячейке F2 тремя знаками после запятой. Вычисленный размах фактических погрешностей составляет 0,007.

Поскольку S вычисляется как интеграл функции одной переменной (размерность интеграла r = 1), то здесь будет справедливо метод Монте-Карло сравнивать с пошаговыми методами на примере метода прямоугольников. Применяя метод прямоугольников (основания которых задаются на отрезках

между точками 0,0000, 0,0001, ... 0,9999), получаем оценку интеграла S ≈ 0,785448. Погрешность этой оценки составляет всего 0,00005. Таким образом, преимущества метода Монте-Карло перед пошаговыми методами начинают проявляться лишь при размерности интеграла r ≥ 2. Заметим, что при размерности r > 2 он становится практически единственным методом, пригодным для численного расчёта интегралов.

В то же время и при размерности r = 1 – в тех случаях, когда область интегрирования имеет сложную конфигурацию, – метод Монте-Карло может оказаться существенно более предпочтительным, чем пошаговые методы.

29

Пример 1 4 . Пусть требуется вычислить определённый интеграл

у которого область интегрирования A определена как совокупность таких непрерывных интервалов на оси t > 0, в которых значения t (в десятичном представлении) имеют сумму s первых трёх цифр после запятой в пределах 3 < s < 18. Например, первый из интервалов области A – это интервал

0,004 ≤ t < 0,010, второй интервал – 0,013 ≤ t < 0,020 и т. д.

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

Выбирая для применения метода в форме (1.46) п. р. в. f(t) = 0,001/t 2, для t > 0,001, мы можем разыгрывать с. в. x, имеющую эту п. р. в., по формуле x = 0,001/z, где z – БСВ [16, с. 42–43]. При попадании с. в. x в область A интегрирования имеем H(x) = h(x) = 1/x2, и, следовательно, усредняемая по методу (1.46) с. в. H(x)/f(x) в области A равна (1/x2) / (1/x2) = 1. Вне области A имеем по определению H(x) = 0. Разыгрывая с. в. x и усредняя H(x), мы получим оценку искомого интеграла I1. Это можно сделать следующим образом.

На чистом листе Ms Excel в диапазоне A1:A10000 сформируйте выборку с. в. x: введите в ячейку A1 формулу =0,001/СЛЧИС() и скопируйте её на весь указанный диапазон ячеек. В соседнем столбце в ячейку B1 введите формулу =A1-ЦЕЛОЕ(A1), вычисляющую дробную часть x, и скопируйте на такой же диапазон ячеек. В следующих трёх столбцах введите в ячейки C1, D1, E1 и аналогично размножьте формулы

=ЦЕЛОЕ(10*B1), =ЦЕЛОЕ(100*B1-10*C1), =ЦЕЛОЕ(1000*B1-100*C1-10*D1),

которые выделяют соответственно первую, вторую и третью цифры после запятой у числа x. Они нужны для того, чтобы определить принадлежность числа x области интегрирования A. В столбце F введите в ячейку F1 формулу вычисле-

ния H(x) (усредняемой с. в.) =ЕСЛИ((C1+D1+E1>3)*(C1+D1+E1<18);1;0) и ско-

пируйте её на диапазон F1:F10000. Определяя далее её среднее и размах погрешностей по правилу трёх сигм, найдите оценку искомого интеграла. Она та-

кова: I1 = 0,22 ± 0,01.

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

4.5. Управление скоростью сходимости

От выбора функции f(t) зависит скорость сходимости оценки интеграла к его точному значению. Скорость сходимости определяется главным образом дисперсией усредняемой в (1.46) с. в. y = H(x)/f(x). Чем меньше D(y), тем мень-

30