12
1.2.2
Сущность метода Монте-Карло
Важнейший прием построения методов Монте-Карло — сведение задачи к расчету математических ожиданий. Пусть требуется найти значение m некоторой изучаемой величины. С этой целью выбирают такую случайную величину X, математическое ожидание которой равно m : M [X] = m.
Практически же поступают так: вычисляют (разыгрывают) N возможных значений x i
случайной величины X, находят их среднее арифметическое
¯
x =
1
N
N
P
i=1
x i
. Так как последовательность одинаково распределенных слу- чайных величин, у которых существуют математические ожидания, под- чиняется закону больших чисел, то при N → ∞ среднее арифметическое этих величин сходится по вероятности к математическому ожиданию. Та- ким образом, при больших N величина ¯
x ≈ m.
В методе Монте-Карло данные вырабатываются искусственно путем ис- пользования некоторого генератора случайных чисел в сочетании с функ- цией распределения вероятностей для исследуемого процесса. Таким ге- нератором может быть таблица, колесо рулетки, подпрограмма ЭВМ или какой-либо другой источник равномерно распределенных случайных чи- сел. Подлежащее разыгрыванию распределение вероятностей может быть основано на эмпирических данных, извлекаемых из ранее сформированных записей, или на результатах последнего эксперимента либо может представ- лять собой известное теоретическое распределение.
Таким образом, для применения метода Монте-Карло необходимо уметь разыгрывать случайную величину.
Приведем так называемые правила для разыгрывания случайных вели- чин.
Введем обозначения: R — непрерывная случайная величина, распреде- ленная равномерно в интервале (0, 1); r i
(i = 1, 2, ...) — случайные числа
(возможные значения R).
Правило 1.
Для того, чтобы разыграть дискретную случайную вели- чину X, заданную законом распределения
X
x
1
x
2
x n
p p
1
p
2
p n
надо:
13
1. Разбить интервал (0, 1) оси Or на n частичных интервалов;
2. Выбрать (например, из таблицы случайных чисел) случайное число r i
Если r i
попало в частичный интервал, то разыгрываемая величина приняла возможное значение x i
Правило 2
(Метод обратных функций). Для того чтобы разыграть возможное значение x i
непрерывной случайной величины X, зная ее функцию распределения F (x), надо выбрать случайное число x i
, прирав- нять его функции распределения и решить относительно x i
полученное уравнение F (x i
) = r i
Правило 3.
Для того, чтобы разыграть возможное значение x i
непре- рывной случайной величины X, зная ее плотность вероятности f (x), надо выбрать случайное число r i
и решить относительно x i
уравнение x
i
Z
−∞
f (x)dx = r i
,
(1.14)
или уравнение x
i
Z
a f (x)dx = r i
(1.15)
где a – наименьшее конечное возможное значение X.
14
1.2.3
Погрешность метода Монте-Карло
Для оценки величины m смоделируем такую случайную величину X,
что ее математическое ожидание M [X] = m. Выберем N независимых реализаций x
1
, . . . , x
N
случайной величины X и вычислим среднее ариф- метическое
¯
x =
1
N
N
X
i=1
x i
Предположим дополнительно, что случайная величина X имеет конеч- ную дисперсию
D[X] = M [X
2
] − (M [X])
2
По центральной предельной теореме:
P
(
| ¯
X − m| < t
β
r
D[X]
N
)
≈ 2Φ(t
β
) = β.
(1.16)
Величина
ε = t
β
r
D[X]
N
(1.17)
называется верхней границей ошибки с коэффициентом доверия β.
Из формулы (1.17) видно, что для того, чтобы уменьшить ошибку в 10
раз (иначе говоря, чтобы получить в ответе еще один верный десятичный знак), нужно увеличить N в 100 раз.
Задаваясь значениями ε и β при известном σ =
pD[X] можно опреде- лить необходимое число испытаний N , обеспечивающих точность ε с на- дежностью β:
N =
σ
2
· t
2
β
ε
2
(1.18)
Однако при решении конкретных задач обычно величина σ оказывается неизвестной. Поэтому определение приближенного значения σ требует до- полнительных рассмотрений. В процессе вычисления можно поступить так.
Назначается ориентировочно предварительное число испытаний N = N
0
Затем по результатам N
0
испытаний определяется приближенное значение дисперсии:
S
2
=
1
N
0
N
0
X
i=1
x
2
i
−
1
N
0
N
0
X
i=1
x i
!
2
(1.19)
15
Располагая S
2
, можно найти приближенное значение N :
N =
S
2
· t
2
β
ε
2
(1.20)
Если окажется, что N > N
0
, необходимо провести дополнительные испы- тания.
16
1.3
Генераторы случайных чисел
Для моделирования случайных процессов, связанных с применением ме- тода Монте-Карло, необходимы случайные числа.
Поскольку при вычислениях методом Монте-Карло существенное коли- чество операций расходуется для оперирования над случайными числами,
то наличие простых и экономных способов формирования последователь- ности случайных чисел во многом определяет возможность практического использования этого метода.
В качестве исходной совокупности случайных чисел, используемых для образования случайных элементов различной природы, необходимо вы- брать такую совокупность, которая может быть получена с наименьши- ми, по возможности, затратами машинного времени и, кроме того, обес- печивает простоту и удобство дальнейших преобразований. Обычно счи- тают, что этим требованиям удовлетворяет совокупность случайных чисел с равномерным распределением в интервале (0, 1). Исходя из равномерно распределенных случайных чисел, можно конструировать как случайные события, возникающие с любой заданной вероятностью, так и случайные величины, обладающие практически любым законом распределения.
Напомним основные свойства равномерного распределения. Непрерыв- ная случайная величина R имеет равномерное распределение в интервале
(a, b), если ее функция плотности равна f (x) =
1
b−a при a ≤ x ≤ b,
0,
вне этого интервала.
(1.21)
Функция распределения случайной величины R имеет вид
F (x) =
0,
при x < a,
x−a b−a при a ≤ x ≤ b,
1,
при x > b.
(1.22)
Математическое ожидание и среднее квадратическое отклонение соответ- ственно равны
M [R] =
a + b
2
,
σ
R
=
b − a
2
√
3
(1.23)
17
В частном случае равномерного распределения в интервале (0, 1):
f (x) =
1 при 0 ≤ x ≤ 1,
0, вне этого интервала,
(1.24)
F (x) =
0, при x < 1,
x при 0 ≤ x ≤ 1,
1, при x > 1,
(1.25)
а математическое ожидание и среднее квадратическое отклонение соответ- ственно равны
M [R] =
1 2
,
σ
R
=
1 2
√
3
(1.26)
Как правило, в качестве стандартной выбирают непрерывную случай- ную величину ξ, равномерно распределенную в интервале (0,1). Заметно реже в качестве стандартной используют дискретную случайную величи- ну η, которая с одинаковой вероятностью может принимать 10 значений
0, 1,..., 9. Мы будем называть величину ξ случайным числом, а величину
η случайной цифрой. Иногда η называют десятичной случайной цифрой,
чтобы отличить ее от двоичной случайной цифры. Чтобы установить связь между ξ и η, разложим число ξ в бесконечную десятичную дробь:
ξ = 0.η
1
η
2
...η
k
Последняя запись означает, что
ξ =
∞
X
k=1
η
k
· 10
−k
Т е о р е м а. Десятичные цифры η
1
, ..., η
k
, ... – случайного числа
ξ представляют собой независимые случайные цифры. Обратно, если
η
1
, ..., η
k
, ... — независимые случайные цифры, то ξ = 0.η
1
η
2
...η
k
... опре- деляет случайное число.
В вычислениях всегда используют числа с конечным количеством деся- тичных знаков, поэтому вместо случайных чисел ξ употребляют конечные десятичные дроби e
ξ = 0.η
1
...η
n
Таблицы случайных цифр. Предположим, что мы осуществили N
независимых опытов, в результате которых получили N случайных цифр
η
1
η
2
...η
N
. Записав эти цифры (в порядке появления) в таблицу, получим то,
что называется таблицей случайных цифр (см. Приложение 1, где цифры
18
объединены в группы только ради удобства чтения). Способ употребления такой таблицы весьма прост. Если в ходе расчета некоторой задачи нам по- требуется случайная цифра η, то мы можем взять любую цифру η
s из этой таблицы. Если нам понадобится случайное число ξ, то мы можем взять из таблицы n очередных цифр и считать, что ξ = 0.η
s
η
s+1
...η
s+n−1
. Выбирать цифры из такой таблицы в случайном порядке не обязательно. Их можно выбирать подряд. Но, конечно, можно начинать с любого места, читать в любом направлении, использовать любой заранее заданный алгоритм вы- бора, не зависящий от конкретных значений цифр таблицы.
Все сказанное выше относится к «идеальной» таблице случайных цифр и не вызывает никаких сомнений. Изготовление хорошей таблицы — весьма сложное дело, ибо в любом реальном опыте всегда возможны ошибки. Про- верка «качества» таблицы абсолютно необходима. Для этого существуют тесты для проверки случайных цифр.
В настоящее время таблицы случайных цифр (или более разнообразные таблицы случайных величин) используют главным образом при расчетах вручную; для расчетов на ЭВМ ими практически не пользуются.
Достоинства метода:
• проверка однократная;
• воспроизводить числа можно.
Недостатки метода:
• запас чисел ограничен;
• занимает много места в накопителе или медленно вводится;
• нужна внешняя память.
Ограниченность объема существующих таблиц, — по-видимому, играет вто- ростепенную роль: можно заготовить таблицу любого объема, если это по- требуется.
Датчики случайных чисел. Чаще всего для построения датчика ис- пользуют «шумящие» радиоэлектронные приборы (диоды, тиратроны, га- зотропы и др.). Не вдаваясь в технические подробности, рассмотрим один из возможных способов построения датчика, вырабатывающего случайные двоичные цифры α. Нетрудно представить себе счетчик, который подсчи- тывает количество ν флуктуаций напряжения шумящего прибора, превы- шающих заданный уровень E
0
за фиксированное время ∆t.
19
s из этой таблицы. Если нам понадобится случайное число ξ, то мы можем взять из таблицы n очередных цифр и считать, что ξ = 0.η
s
η
s+1
...η
s+n−1
. Выбирать цифры из такой таблицы в случайном порядке не обязательно. Их можно выбирать подряд. Но, конечно, можно начинать с любого места, читать в любом направлении, использовать любой заранее заданный алгоритм вы- бора, не зависящий от конкретных значений цифр таблицы.
Все сказанное выше относится к «идеальной» таблице случайных цифр и не вызывает никаких сомнений. Изготовление хорошей таблицы — весьма сложное дело, ибо в любом реальном опыте всегда возможны ошибки. Про- верка «качества» таблицы абсолютно необходима. Для этого существуют тесты для проверки случайных цифр.
В настоящее время таблицы случайных цифр (или более разнообразные таблицы случайных величин) используют главным образом при расчетах вручную; для расчетов на ЭВМ ими практически не пользуются.
Достоинства метода:
• проверка однократная;
• воспроизводить числа можно.
Недостатки метода:
• запас чисел ограничен;
• занимает много места в накопителе или медленно вводится;
• нужна внешняя память.
Ограниченность объема существующих таблиц, — по-видимому, играет вто- ростепенную роль: можно заготовить таблицу любого объема, если это по- требуется.
Датчики случайных чисел. Чаще всего для построения датчика ис- пользуют «шумящие» радиоэлектронные приборы (диоды, тиратроны, га- зотропы и др.). Не вдаваясь в технические подробности, рассмотрим один из возможных способов построения датчика, вырабатывающего случайные двоичные цифры α. Нетрудно представить себе счетчик, который подсчи- тывает количество ν флуктуаций напряжения шумящего прибора, превы- шающих заданный уровень E
0
за фиксированное время ∆t.
19
Рис. 1.1: Флуктуации напряжения шумящего прибора
Еще проще устроить счетчик, который выдавал бы число ν(mod2), т.
е. 0 при четном ν и 1 при нечетном ν. Если вероятности появления 0 и 1
в таком процессе равны между собой, то можно считать, что устройство вырабатывает случайную последовательность двоичных цифр.
Обычно датчики случайных чисел содержат m генераторов описанного типа, работающих независимо, так что датчиком выдается приближенное случайное число ξ = 0.α
1
...α
m
, записанное в форме m-разрядной двоичной дроби. Для случайных чисел отведена специальная ячейка в накопителе,
и скорость генерирования их столь велика, что на каждом такте работы
ЭВМ в этой ячейке получается новое случайное число.
Достоинства метода:
• запас чисел неограничен;
• сверхбыстрое получение;
• места в накопителе не занимает.
Недостатки метода:
• проверка периодическая;
• воспроизводить числа нельзя;
• требуется специальное устройство.
Псевдослучайные числа. Пригодность случайных чисел определя- ется в конечном счете не процессом их получения, а тем, удовлетворяют ли они некоторым принятым тестам. Но в таком случае совершенно без- различно, как эти числа получены, они могут быть даже сосчитаны по какой-нибудь формуле, лишь бы они удовлетворяли тестам.
20
Числа ξ
1
, ξ
2
, ..., ξ
N
, которые вычисляются по какой-либо заданной фор- муле и могут быть использованы вместо случайных чисел при решении некоторых задач, называются псевдослучайными числами.
Большинство алгоритмов, используемых на практике для получения псевдослучайных чисел, представляют собой рекуррентные формулы пер- вого порядка
ξ
n+1
= F (ξ
n
),
(1.27)
где начальное число ξ
0
задано.
Важной чертой алгоритмов вида (1.27) является то, что при реализации их на ЭВМ они всегда порождают периодические последовательности.
Достоинства метода:
• проверка однократная;
• воспроизводить числа можно;
• быстрое получение;
• места в накопителе занимает мало;
• внешние устройства не нужны.
Недостаток – запас чисел ограничен.
Метод псевдослучайных чисел – самый удобный с практической точ- ки зрения. Это подтверждается также всей практикой расчетов методами
Монте-Карло.
21
1.4
Интегрирование с помощью метода Монте-Карло
1.4.1
Способ усреднения подынтегральной функции
Рассмотрим интеграл вида
I =
b
Z
a
ϕ(x)dx.
(1.28)
Введем в рассмотрение случайную величину X, распределенную равно- мерно в интервале интегрирования (a, b) с плотностью f (x) =
1
b−a
. Тогда математическое ожидание
M [ϕ(x)] =
b
Z
a
ϕ(x)f (x)dx =
1
b − a
Z
b a
ϕ(x)dx.
(1.29)
Отсюда b
Z
a
ϕ(x)dx = (b − a)M [ϕ(x)].
(1.30)
Заменив математическое ожидание M [ϕ(x)] его оценкой – выборочной средней, получим оценку интеграла (1.28)
I
∗
1
= (b − a)
N
P
i=1
ϕ(x i
)
N
,
(1.31)
где x i
– возможные значения случайной величины X, N – число испытаний.
Так как случайная величина распределена равномерно в интервале (a, b) с плотностью f (x) =
1
b−a
, то x i
разыгрываются по формуле
1
b − a x
i
Z
a dx = r i
(1.32)
Отсюда x i
= a + (b − a)r i
, где r i
- случайное число.
22
Таким образом, в качестве оценки определенного интеграла (1.28) при- нимают (1.31).
Дисперсия усредняемой функции ϕ(X) равна
σ
2
=
b
Z
a
ϕ
2
(x)f (x)dx − (M [ϕ(X)])
2
(1.33)
Если точное значение дисперсии вычислить трудно или невозможно, то находят выборочную дисперсию (при N > 30),
S =
N
X
i=1
u
2
i
N
−
N
X
i=1
u i
N
!
2
,
(1.34)
или исправленную дисперсию (при N < 30)
S =
1
N − 1
·
N
X
i=1
u
2
i
N
−
N
P
i=1
u i
2
N
,
(1.35)
где u i
= ϕ(x i
).
Эти формулы для вычисления дисперсии применяют и при других способах интегрирования, когда усредняемая функция не совпадает с подынтегральной функцией.
Пример 1.1
Найти:
а) оценку определенного интеграла I =
3
R
1
(x + 1)dx;
б) абсолютную погрешность |I − I
∗
|;
в) минимальное число испытаний, которые с надежностью 0.95 обеспе- чат верхнюю границу ошибки 0.1.
Решение.
а) Используем формулу
I
∗
1
= (b − a)
N
P
i=1
ϕ(x i
)
N
23
По условию, a = 1, b = 3, ϕ(x) = x + 1. Примем для простоты число испытаний N = 10. Тогда оценка I
∗
1
= 2 ·
10
P
i=1
(x i
+1)
10
, где возможные значения x
i разыгрываются по формуле x
i
= a + (b − a) · r i
= 1 + 2 · r i
Результаты десяти испытаний приведены в таблице. Случайные числа r
i взяты из таблицы Приложения 1.
Номер испытания, i
1 2
3 4
5 6
7 8
9 10
r i
0.100 0.973 0.253 0.376 0.520 0.135 0.863 0.467 0.354 0.876
i + 2r i
1.200 2.946 1.506 1.752 2.040 1.270 2.726 1.934 1.708 2.752
x i
+ 1 2.200 3.946 2.506 2.752 3.040 2.270 3.726 2.934 2.708 3.752
Из таблицы находим
10
P
i=1
ϕ(x i
) = 29.834. Искомая оценка
I
∗
1
= 2 ·
29.834 10
= 5.967
б) Приняв во внимание, что I =
3
R
1
(x + 1)dx = 6, найдем абсолютную погрешность
|6 − 5.967| = 0.033
в) Найдем дисперсию усредняемой функции ϕ(X) = X + 1, учитывая,
что случайная величина X в интервале интегрирования (1, 3) распределена равномерно и ее дисперсия
D[X] = (3 − 1) ·
2 12
=
1 3
σ
2
= D[X + 1] = D[X] =
1 3
Найдем минимальное число испытаний, которые с надежностью 0.95
обеспечат верхнюю границу ошибки ε = 0, 1. Из равенства Φ(t) = 0.95/2 =
0.475 по таблице Приложения 2 находим t = 1.96. Искомое минимальное число испытаний
N =
t
2
· σ
2
ε
2
=
1.96 2
·
1 3
0.1 2
= 128.
24
Пример 1.2
В качестве приближенного значения интеграла I =
π/2
R
0
cos xdx принята оценка (1.31). Найти минимальное число испытаний, при котором с надежностью 0.95 верхняя граница ошибки δ = 0.1.
Решение.
Найдем дисперсию усредняемой функции ϕ(x) = cos X,
учитывая, что случайная величина X распределена равномерно в интер- вале (0,
π
2
) с плотностью f (x) =
π
2
:
σ
2
=
π
2
π/2
Z
0
cos
2
xdx −
π
2
π/2
Z
0
cos xdx
2
= 0.09.
Из равенства Φ(t) = 0.95/2 = 0.475 по таблице приложения (2.7) нахо- дим t = 1.96. Искомое минимальное число испытаний
N =
t
2
· σ
2
ε
2
=
1, 96 2
· 0, 09 0, 1 2
= 35.
25
1.4.2
Способ существенной выборки, использующей вспомога- тельную плотность распределения
Для вычисления интеграла (1.28) введем функцию f (x) - плотность рас- пределения вспомогательной случайной величины X в интервале интегри- рования (a, b), т.е.
b
Z
a f (x)dx = 1.
(1.36)
Представим (1.28) так
I =
b
Z
a
ϕ(x)
f (x)
· f (x)dx.
(1.37)
I представлен в виде математического ожидания функции F (X) =
ϕ(X)
f (X)
В качестве оценки этого математического, а следовательно равного ему интеграла (1.28), примем выборочную среднюю
I
∗
2
=
1
N
N
X
i=1
ϕ(x i
)
f (x i
)
(1.38)
Отметим, что желательно выбрать f (x) так, чтобы по возможности от- ношение
ϕ(x)
f (x)
было постоянным, поскольку в этом случае реализуется ми- нимальная дисперсия, т.е. погрешность будет минимальна ([3], стр. 109.).
Таким образом, в качестве оценки интеграла (1.28) принимают (1.38),
где f (x) – плотность распределения вспомогательной случайной величины
X, где x i
– возможные значения случайной величины X, которые разыг- рываются по формуле x
i
R
a f (x)dx = r i
; N – число испытаний.
Пример 1.3
1 2 3 4 5