Лоп^лярныс лекции
ПО МАТЕМАТИКЕ
И.М.СОБОЛЬ
МЕТОД
МОНТЕ-КАРЛО

книги
МКИТИНА

ПОПУЛЯРНЫЕ ЛЕКЦИИ ПО МАТЕМАТИКЕ
ВЫПУСК 46
И. М. СОБОЛЬ
МЕТОД
МОНТЕ-КАРЛО
ИЗДАТЕЛЬСТВО «НАУКА»
ГЛАВНАЯ РЕДАКЦИЯ
ФИЗИКО-МАТЕМАТИЧЕСКОЙ ЛИТЕРАТУРЫ
МОСКВА 1968

517.8 + 518
С 54
УДК 519.2
Соболь И. М.
С 54	Метод Монте-Карло. М., «Наука», 1968,
64 с. («Популярные лекции по математике», вып. 46), 10 к.
В книжке изложены основные приемы метода Монте-Карло
(метода статистических испытаний). Приведены примеры весьма
разнообразных задач, решаемых этим методом. Предназначена
для инженеров, конструкторов, исследователей и научных работ-
работников, работающих в различных отраслях народного хозяйства
(в науке, технике, промышленности, медицине, экономике, сель-
сельском хозяйстве, торговле и др.), н ставит своей целью подска-
подсказать, что в любой области деятельности встречаются задачи, ко-
которые можно рассчитывать методом Монте-Карло. К читателю
предъявляются минимальные требования: умение дифференциро-
дифференцировать и интегрировать A-й курс втуза). Книжка может быть по-
полезна также всем, кто желает впервые познакомиться с методом
Монте-Карло.
2-2-2	517.8 + 518
147-68
Илья Меерович Соболь
Метод Монте-Карло
М., 1968 г., 64 стр. с илл.
Редакторы М. М. Горячая и И. М. Овчинникова
Техн. редактор В. Н. Крючкова	Корректор Н. В, Гераськина
Сдано в набор 24/Х 1967 г. Подписано к печати 21/Н 1968 г. Бумага 84х108'/м.
Физ. печ. л. 2. Условн. печ. л. 3,36. Уч.-изд. л. 2,90. Тираж 79 000 экз. Т-00260.
Цена 10 коп. Заказ № 925.
Издательство «Наука»
Главная редакция Физико-математической литературы
Москва, В-71, Ленинский проспект, 15.
Ленинградская типография № 2 имени Евгении Соколовой
Главполнграфпрома Комитета по печати при Совете Министров СССР.
Измайловский проспект, 29,

ПРЕДИСЛОВИЕ
Несколько лет назад мне было предложено прочесть
две лекции о методе Монте-Карло для слушателей фа-
факультета вычислительной техники Общественного уни-
университета. Я согласился. И перед первой лекцией с ужа-
ужасом узнал, что большая часть слушателей не знакома с
теорией вероятностей. Отступать было поздно. Пришлось
на ходу включить в лекцию дополнительный раздел, зна-
комящий слушателей с основными теоретико-вероятност-
теоретико-вероятностными понятиями.
Лекции эти повторялись в течение нескольких лет.
Постепенно содержание их «устоялось». В настоящем
издании сохранен также дополнительный раздел (§ 2),
о котором хочется сказать несколько слов.
Каждому человеку приходилось употреблять слова
«вероятность» и «случайная величина». Интуитивное
представление о вероятности (как о частоте) более или
менее соответствует истинному смыслу этого понятия.
Но интуитивное представление о случайной величине,
как правило, весьма далеко от математического опреде-
определения. Поэтому в § 2 понятие вероятности предполагает-
предполагается известным и разъясняется только более сложное по-
понятие случайной величины. Этот параграф не может за-
заменить курса теории вероятностей: изложение здесь
упрощенное и без доказательств. Но он дает некоторое
представление о случайных величинах, достаточное для
понимания простейших приемов метода Монте-Карло.
1*	3

Основная цель настоящей книжки — подсказать спе-
специалистам самых различных направлений, что в их об-
области деятельности встречаются задачи, которые можно
решать методом Монте-Карло.
Задачи, рассмотренные в лекциях, достаточно просты
и разнообразны. Но они, конечно, не могли охватить
всех областей применения метода. Ограничусь только
одним примером. В книжке ни слова нет о медицине.
Однако методы § 7 позволяют рассчитывать дозы облу-
облучения при лучевой терапии. Имея программу для рас-
расчета поглощенного излучения в различных тканях тела,
можно наиболее эффективным образом выбирать дози-
дозировку и направление облучения, следя за тем, чтобы не
нанести вреда здоровым тканям.
В настоящую книжку включено все то, что читалось
в лекциях, несколько подробнее изложены примеры и
добавлен § 9.
И. Соболь
Москва, 1967 г.

ВВЕДЕНИЕ
§ 1. Общее представление о методе
Метод Монте-Карло — это численный метод решения
математических задач при помощи моделирования слу-
случайных величин.
1.1. Происхождение метода Монте-Кар-
л о. Датой рождения метода Монте-Карло принято счи-
считать 1949 год, когда появилась статья под названием
«The Monte Carlo method»*). Создателями этого метода
считают американских математиков Дж. Неймана и
С. Улама. В Советском Союзе первые статьи о методе
Монте-Карло были опубликованы в 1955—1956 годах**)'.
Любопытно, что теоретическая основа метода была
известна уже давно. Более того, некоторые задачи ста-
статистики рассчитывались иногда с помощью случайных
выборок, то есть фактически методом Монте-Карло.
Однако до появления электронных вычислительных ма-
машин (ЭВМ) этот метод не мог найти сколько-нибудь ши*
рокого применения, ибо моделировать случайные вели*
чины вручную —очень трудоемкая работа. Таким обра»
зом, возникновение метода Монте-Карло как весьма
универсального численного метода стало возможным
только благодаря появлению ЭВМ.
Само название «Монте-Карло» происходит от города
Монте-Карло в княжестве Монако, знаменитого своим
игорным домом. Дело в том, что одним из простейших
механических приборов для получения случайных величин
*) Metropolis N., Ulam S., The Monte Carlo method, J.
Amer. statistical assoc, 1949, 44, № 247, 335—341.
**) Это были статьи В. В. Чавчанидзе, Ю. А. Шрейдера,
В. С. Владимирова.

является... рулетка. Об этом речь пойдет в § 3. Здесь
же, пожалуй, стоит ответить на часто задаваемый во-
вопрос: «Помогает ли метод Монте-Карло выигрывать в
рулетку?» Нет, не помогает. И даже не занимается этим.
1.2.	Пример. Для того чтобы читателю было более
понятно, о чем пойдет речь, рассмотрим весьма простой
пример. Предположим, что нам нужно вычислить пло-
площадь плоской фигуры S. Это
может быть совсем произ-
произвольная фигура с криволи-
криволинейной границей, заданная
графически или аналитиче-
аналитически, связная или состоящая
из нескольких кусков. Пусть
это будет фигура, изобра-
изображенная на рис. 1, и предпо-
предположим, что она вся распо-
расположена внутри единичного
квадрата.
Выберем в квадрате N
х случайных точек. Обозначим
Рис. 1.	через N' число точек, попав-
попавших при этом внутрь S. Гео-
Геометрически очевидно, что площадь 5 приближенно равна
отношению N'/N. Чем больше будет iV, тем больше будет
точность этой оценки.
В примере, изображенном на рис. 1, выбраны jV = 40
точек. Из них N'=\2 точек оказались внутри S. Отноше-
Отношение N'IN= 12/40 = 0,30, в то время как истинная пло-
площадь S равна 0,35 *).
1.3.	Две особенности метода Монте-Кар-
Монте-Карло. Первая особенность метода — простая структура вы-
вычислительного алгоритма. Как правило, составляется
программа для осуществления одного случайного испы-
испытания (в примере из п. 1.2 надо выбрать случайную
*) На практике для вычисления площади плоской фигуры метод
Монте-Карло не используют: для этого есть другие методы, хотя и
более сложные, но зато обеспечивающие гораздо большую точность.
Однако указанный в нашем примере метод Монте-Карло позво-
позволяет столь же просто вычислять «многомерный объем» тела в мно-
многомерном пространстве. И в этом случае метод Монте-Карло часто
оказывается единственным численным методом, дающим возмож-
возможность решнть задачу,
б

точку в квадрате и проверить, принадлежит ли она 5).
Затем это испытание повторяется N раз, причем каждый
опыт не зависит от всех остальных, и результаты всех
опытов осредняются.
Поэтому иногда метод Монте-Карло называют мето*
дом статистических испытаний.
Вторая особенность метода: ошибка вычислений, как
правило, пропорциональна У DIN, где D — некоторая по-
постоянная, а N — число испытаний. Из этой формулы вид*
но, что для того, чтобы уменьшить ошибку в 10 раз
(иначе говоря, чтобы получить в ответе еще один вер-
ный десятичный знак), нужно увеличить N (то есть
объем работы) в 100 раз.
Ясно, что добиться высокой точности на таком пути
невозможно. Поэтому обычно говорят, что метод Монте-
Карло особенно эффективен при решении тех задач, в
которых результат нужен с небольшой точностью
,E-10%).
Однако одну и ту же задачу можно решать различи
ными вариантами метода Монте-Карло*), которым отве-
отвечают различные значения D. Во многих задачах удается
значительно увеличить точность, выбрав способ расчета,
которому соответствует значительно меньшее значение/).
1.4. Задачи, решаемые методом Монте-
Карло. Во-первых, метод Монте-Карло позволяет моде-
моделировать любой процесс, на протекание которого влияют
случайные факторы. Во-вторых, для многих математиче-
математических задач, не связанных с какими-либо случайностями,
можно искусственно придумать вероятностную модель
(и даже не одну), позволяющую решать эти задачи. Соб-
Собственно говоря, это и было сделано в примере из п. 1.2.
Таким образом, можно говорить о методе Монте-
Карло как об универсальном методе решения математи-
математических задач.
Особенно интересно, что в некоторых случаях выгод-
выгодно отказаться от моделирования истинного случайного
процесса и вместо этого использовать искусственную
модель. Такая ситуация рассмотрена в § 7.
*) В иностранной литературе теперь чаще пишут о методах
Монте-Карло (во множественном числе), имея в виду, что одну и
ту же задачу можно считать при помощи моделирования различных
случайных величин.

1.5. Еще о п р и м ер е. Вернемся к примеру из п. 1.2.
Для расчета нам нужно было выбирать случайные точки
в единичном квадрате. Как это сделать физически?
Представим такой эксперимент. Рис. 1 (в увеличен-
увеличенном масштабе) с фигурой S и квадратом повешен на
стену в качестве мишени. Стрелок, находящийся на не-
некотором расстоянии от стены, стреляет N раз, целясь в
центр квадрата. Конечно, все пули не будут ложиться
точно в центр: они про-
пробьют на мишени N слу-
случайных точек*). Можно
ли по этим точкам оце-
оценить площадь S?
Результат такого опы-
опыта изображен на рис. 2.
В этом опыте N = 40, N'=*
= 24 и отношение N'IN =
= 0,60, что превышает
истинное значение площа-
площади @,35) почти в два раза.
Впрочем, и так ясно, что
при очень высокой квали-
Рнс. 2.	фикации стрелка резуль-
результат опыта будет очень
рлохим, так как почти все пули будут ложиться вблизи
центра и попадут в S**).
Нетрудно сообразить, что наш метод вычисления пло-
площади будет справедлив только тогда, когда случайные
точки будут не «просто случайными», а еще и «равно-
«равномерно разбросанными» по всему квадрату. Чтобы при-
придать этим словам точный смысл, необходимо познако-
познакомиться с определением случайных величин и некоторыми
их свойствами. Эти сведения изложены в § 2 без дока-
доказательств. Читатель, знакомый с любым курсом теории
вероятностей, может пропустить весь § 2, кроме п. 2.5.
*) Мы предполагаем, что стрелок не является чемпионом мира
и удален на достаточно большое расстояние от мишени,
**) О том, как выбирались случайные точки на рис, 1 и 2,
будет сказано в п. 4.5.

ГЛАВА 1
МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ВЕЛИЧИН
§ 2. Случайные величины
Мы считаем, что понятие вероятность читателю более
или менее известно, и переходим прямо к понятию слу«
чайная величина.
Слова «случайная величина» в обыденном смысле
употребляют тогда, когда хотят подчеркнуть, что неиз-
неизвестно, каким будет конкретное значение этой величины,
Причем иногда за этими словами скрывается просто не«
знание, какова эта величина.
Однако математик употребляет эти же слова «слу»
чайная величина», вкладывая в них вполне определен*
ный положительный смысл. Действительно, говорит он,
мы не знаем, какое значение примет эта величина в дан-
данном конкретном случае, но мы знаем, какие значения
она может принимать, и знаем, каковы вероятности тех
или иных значений. На основании этих данных мы не
мож«м точно предсказать результат одного испытания,
связанного с этой случайной величиной, но можем весь-
весьма надежно предсказать совокупность результатов боль-
большого числа испытаний. Чем больше испытаний (как го-
говорят, чем больше статистика), тем точнее будут наши
предсказания.
Итак, чтобы задать случайную величину, надо ука*
зать, какие значения она может принимать и каковы
вероятности этих значений.
2.1. Дискретные случайные величины. Слу«
чайная величина \ называется дискретной, если она
9

может принимать дискретное множество значений xit
Х2, • • •> Хп ) .
Дискретная случайная величина | определяется таб-
таблицей
Лу А2 . . . Лл
(Т)
i A ••• /»я/
где хь а'2, ..., а'„ — возможные значения величины |, а
/?ь /?2, ..., Рп — соответствующие им вероятности. Точ-
Точнее говоря, вероятность того, что случайная величина |
примет значение Х\ (обозначим ее через Pi^—Xi)), рав-
равна П:\
Таблица (Т) называется распределением случайной
величины.
Числа Х\, х2, ..., хп могут быть, вообще говоря, лю-
любыми. Однако вероятности р1} рг, ¦.., рп должны удо-
удовлетворять двум условиям:
а)	все pi положительны:
Pi > 0;	A)
б)	сумма всех pt равна 1:
А + А+ ••¦ +/>„=!•	B)
Последнее условие означает, что \ обязана в ка-
каждом случае принять одно из значений xit х2, ..., хп.
Математическим ожиданием случайной величины |
называется число
щ=i xiPi.	C)
Чтобы выяснить физический смысл этой величины,
запишем ее в следующем виде:
mi-
*) В теории вероятностей рассматриваются также дискретные
случайные величины, которые могут принимать бесконечное число
значений.
10

Отсюда видно, что М| — это среднее значение величи-
величины |, причем более вероятные значения хг входят в сум-
сумму с большими весами*).
Отметим основные свойства математического ожида*
ния: если с — какая-нибудь не случайная величина, то
D)
М (с\) = сМ|;	E)
если | и г) — две любые случайные величины, то
M(S + ti) = M|+Mti.	F)
Дисперсией случайной величины | называется	число
Di = M[(i-M?J].	G)
Следовательно, дисперсия D1 — это математическое ожи->
дание квадрата отклонения случайной величины | от ее
среднего значения М|. Очевидно, всегда DE> 0.
Математическое ожидание и дисперсия — важнейшие
числовые характеристики случайной величины |. Каково
их практическое значение?
Если мы будем наблюдать величину | много раз и
получим значения |4, |2, • •., In (каждое из которых бу-
будет равно одному из чисел хи х2, ..., хп), то среднее
арифметическое от этих значений будет близко к М|
¦^(Ii + E2+ ... +Ы^М|.	(8)
А дисперсия D? характеризует разброс этих значений
около среднего М|.
*) Осреднение с весами очень часто встречается в самых раз-
различных областях науки. Например, в механике: если в точках
*ь *2, ..., Хп (на оси Ох) расположены массы mb mi, ..,, тп, то
абсцисса центра тяжести этой системы вычисляется по формуле
(-1
Конечно, в этом случае сумма всех масс не обязана равняться еди-
единице.
11

Формулу G) для дисперсии можно преобразовать е
помощью формул D) — F):
D1 = М [|2 - 2MS ¦ | + (М|J] = М (|2) - 2MS • М| + (MEJ,
откуда
DE = M(E2)-(M|J.	(9)
Вычислять дисперсию по формуле (9) обычно проще,
чем по формуле G).
Отметим основные свойства дисперсии: если с —
какая-нибудь не случайная величина, то
D(S+c) = DS,	(Ю)
D(cl) = c°Dl.	(И)
Важную роль в теории вероятностей играет понятие
независимости случайных величин. В действитель-
действительности это — довольно сложное понятие, но в простейших
случаях оно весьма очевидно*).
Для независимых случайных величин | и г| справед-
справедливы следующие соотношения:
М(?тО = М|.Мть	A2)
A3)
Пример. Рассмотрим случайную величину к с распределе-
распределением
/1 2 3 4 5 6 \
Х \1/6 1/6 1/6 1/6 1/6 1/6/
Очевидно, реализацией этой величины можно считать число очков,
выпадающих на игральной кости: любое значение одинаково веро-
вероятно. Вычислим математическое ожидание и дисперсию ч. По фор-
формуле C)
Му. = ~-A+2 + 3 + 4 + 5+6) =3,5.
По формуле (9)
Dx = М (к2) — (МхJ = -1 (I2-f-22 + ... +62) — C,5J = 2,917.
*) Допустим, что кроме величины | мы наблюдаем еще за од-
одной случайной величиной т]. Если распределение величины | не ме-
меняется от того, что нам уже известно значение, которое приняла
величина т\, то естественно считать, что ^ от т] не зависит,
12

Пример. Рассмотрим случайную величину 9 с распределе-
нием
/2 1/2/
Реализацией этой величины можно считать игру в орлянку с
условием, что, например, орел — это 3 очка, а решетка — 4 очка.
В stom случае
МО = 0,5 • 3 -f- 0,5 • 4 = 3,5;
DG = 0,5 (З2 + 42) — C.5J = 0,25.
Мы видим, что MG = Мх, но DO < Dx. Это можно было легко
предвидеть, так как значения 6 могут отличаться от 3,5 только иа
±0,5, в то время как в значениях х разброс может достигать ±2,5.
2.2. Непрерывные случайные величины.
Предположим, что иа плоскости в начале координат рас-
расположено некоторое количество радия. При распаде
каждого атома радия из него вылетает а-частица. На-
Направление ее будем характери-
характеризовать углом 1)з (рис. 3). Так как ,,
и теоретически, и практически
возможны любые направления
вылета, то эта случайная вели-
величина может принимать любое
значение от 0 до 2я.
Мы будем называть случай-
случайную величину I непрерывной, если _
она может принимать любое зна- о
чение из некоторого интервала
(а, Ь).	Рис. 3.
Непрерывная случайная ве-
величина | определяется заданием
интервала (а, Ь), содержащего возможные значения
этой величины, и функции р(х), которая называется
плотностью вероятностей случайной величины | (или
плотностью распределения g).
Физический смысл р(х) следующий: пусть (а', Ь') —
произвольный интервал, содержащийся в (а, Ь) (то есть
а 4^. a', b'-4ib). Тогда вероятность того, что | окажется
в интервале (а', Ь'), равна интегралу
Р[а'
= \ p{x)dx.
A4)

На рис. 4 заштрихованная площадь равна значению ин-
интеграла A4).
Множество значений | может быть любым интерва-
интервалом. Возможен даже случай а ——оо, а также Ь = оо.
a' b	b л
Рис. 4.
Однако плотность р(х) должна удовлетворять двум ус-
условиям, аналогичным условиям A) и B) для дискрет-
дискретных величин:
а)	плотность р(х) положительна:
/>(*)> 0;	A5)
б)	интеграл от плотности р(х) по всему интервалу
(а, Ь) равен 1:
ь
p(x)dx=l.	A6)
Математическим ожиданием непрерывной случайной
величины называется число
ь
щ = J xp (x) dx.	A7)
а
Смысл этой характеристики такой же, как в случае
дискретной случайной величины. В самом деле, так как
J хр (х) dx
i dx
14

то легко видеть, что это есть среднее значение |: ведь
значением | может быть любое число х из интервала
(й, Ь), которое входит в интеграл с весом р(х)*).
Все изложенное в п. 2.1 от формулы D) и до фор-
формулы A3) включительно справедливо также для непре*
рывных случайных величин: и определение дисперсии
G), и формула (9) для ее вычисления, и все свойства
М| и D|- Повторять мы их не будем.
Укажем только еще одну формулу для математиче*
ского ожидания случайной функции. Пусть по-прежнему
случайная величина | имеет плотность вероятностей
р(х). Выберем произвольную непрерывную функцию f(x)
и рассмотрим случайную величи-
величину "Л —/(?). которую иногда на-
зывают случайной функцией.
Можно доказать, что
ь
Mf(l)=lf(x)p(x)dx. A8)
Подчеркнем, что, вообще гово- л
ря, Щ)П
Случайная величина у, опре-
определенная в интервале @, 1) и
имеющая плотность р(х) = \, называется равномерно
распределенной в (О, 1) (рис. 5).
В самом деле, какой бы интервал (й', Ь') внутри
(О, 1) мы ни взяли, вероятность того, что у попадет в
(V, Ь'), равна
Ь'
J p{x)dx = b' —а',
а'
*) И в этом случае можно указать на аналогичную формулу
из механики: если линейная плотность стержня а< л-< b равна
р(х), то абсцисса центра тяжести стержня вычисляется по формуле
J хр (х)
dx
— а
\dx
15

то есть длине этого интервала. В частности, если мы
разделим @, 1) на любое число интервалов равной дли*
ны, то вероятность попадания у в любой из этих интер-
интервалов будет одна и та же.
Легко вычислить, что
My == f хр (х) dx = Г х dx — ^;
о	о
Со случайной величиной у мы еще не раз встретимся
в дальнейшем.
2.3. Нормальные случайные величины.
Нормальной (или гауссовской) случайной величиной на-
называется случайная величина ?, определенная на всей
оси (—с», оо) и имеющая плотность
где а и о>0 — числовые параметры*).
Параметр а не влияет на форму кривой р(х): изме-
изменение его приводит лишь к сдвигу кривой вдоль оси х.
Однако при изменении о форма кривой меняется. В са-
мом деле, легко видеть, что
max р(х) = р (а) = ——=-.
Если уменьшать о, то тахр(лг) будет возрастать. Однако
вся площадь под кривой р(х) по условию A6) равна 1.
Поэтому кривая будет вытягиваться вверх в окрестности
х = а, но убывать при всех достаточно больших значе-
значениях х. На рис. 6 построены две нормальные плотности,
соответствующие а=0, о=1 и я = 0, о=0,5. (Еще одна
нормальная плотность построена на рис. 21, см. ниже,
стр. 43.)
Можно доказать, что
*) а — это число, а не случайная величина. Мы использовали
греческую букву, так как это обозначение общепринято.
16

Нормальные случайные величины очень часто встре-
встречаются при исследовании самых различных по своей при-
природе вопросов. О причине этого будет сказано несколько
ниже. Например, ошибка измерения б, как правило,
представляет собой нормальную случайную величину.
Если систематической
=М6 — 0. А величина о
квадратичной погрешностью,
ошибки измерения нет, то й=?
= YD&, называемая средней
характеризует
погрешность метода измерения.
Правило «трех сиг м». Нетрудно вычислить, что,
каковы бы ни были аи о в A9),
в-:3о
f n(r\/1\	О QQ7
Из A4) вытекает, что
Р {а — Зо < I < a -f Зо} = 0,997.
B0)
Вероятность 0,997 настолько близка к 1, что иногда
последнюю формулу интерпретируют так: при одном
испытании практически невозможно по-
получить значение ?, отличающееся от MS
больше чем на Зо.
2.4. Центральная предельная теорема
теории вероятностей. Эта замечательная теорема
2 И. М. Соболь
17

была впервые сформулирована Лапласом. Обобщением
этой теоремы занимались многие выдающиеся матема-
математики, в том числе П. Л. Чебышев, А. А. Марков и А. М. Ля-
Ляпунов. Доказательство ее достаточно сложно.
Рассмотрим N одинаковых независимых случайных
величин |ь |2, .. •. liv. Иначе говоря, плотности вероят-
вероятностей этих величин совпадают. Следовательно, их мате-
математические ожидания и дисперсии также совпадают.
Обозначим
Обозначим через рл* сумму всех этих величин:
PJv==:bi + E2 4- ••• -{-In-
Из формул F) и A3) следует, что
Nip* = М (|, + S2 4-	1- In) = Nm>
Рассмотрим теперь нормальную случайную величину
t,N с такими же параметрами: a~Nm, o2 = Nb2.
Теорема. Плотность суммы pN приближается к
плотности нормальной величины t,N, когда N -> оо:
Физический смысл этой теоремы очевиден: сумма pv
большого числа одинаковых случайных величин прибли-
приблизительно нормальна (рр (х) ^ pt (x) \.
На самом деле эта теорема справедлива при гораздо
более широких условиях: все слагаемые |ь |2, ..., |,v
не обязаны быть одинаковыми и независимыми; суще-
существенно только, чтобы отдельные слагаемые не играли
слишком большой роли в сумме.
Именно эта теорема объясняет, почему нормальные
случайные величины так часто встречаются в природе.
В самом деле, каждый раз, когда мы сталкиваемся с
Суммарным воздействием большого числа незначитель-
незначительных случайных факторов, результирующая случайная
величина оказывается нормальной.
18

Например, отклонение артиллерийского снаряда от
цели почти всегда оказывается нормальной случайной
величиной, так как оно зависит и от метеорологических
условий на разных участках траектории, и от многих
других факторов.
2.5. Общая схема метода Монте-Карло.
Допустим, что нам требуется вычислить какую-то неиз-
неизвестную величину т. Попытаемся придумать такую слу-
случайную величину |, чтобы Ы% — т. Пусть при этом
Рассмотрим N случайных величин |i, |2, • ¦., \n, рас-
распределения которых совпадают с распределением |.
Если N достаточно велико, то согласно теореме из п. 2.4
распределение суммы pjv = 5i + l2+. ¦ . + |jv будет прибли-
приблизительно нормальным с параметрами a = Nm, G***Nb\
Из B0) следует, что
Р [Nm — ЗЬ УЫ < pjV < Nm + 3b ]/7v} =^0,997.
Если мы поделим неравенство, стоящее в фигурной
скобке, на N, то получим эквивалентное неравенство и
вероятность его останется такой же:
Pirn	-p=r < .!&_ < m + l я** o,997.
Последнее соотношение перепишем в несколько ином
виде:
1
¦т
3ft
Это — чрезвычайно важное для метода Монте-Карло
соотношение. Оно дает нам и метод расчета т, и оценку
погрешности.
В самом деле, найдем N значений случайной величи-
величины |*). Из B1) видно, что среднее арифметическое
*) Все равно, находить ли один раз по одному значению каж-
каждой из величин |i, |г	|n иди найти iV значений одной вели-
величины |, так как все эти случайные величины совпадакй (имеют
одно и то же распределение),
2*	19

этих значений будет приближенно равно т. С большой
вероятностью ошибка такого приближения не превосхо-
превосходит величины ZbJYN. Очевидно, эта ошибка стремится
к нулю с ростом N.
§ 3. Получение случайных величин на ЭВМ
Сама постановка вопроса «получение случайных ве-
величин на ЭВМ» иногда вызывает недоумение: ведь все,
что делает машина, должно быть заранее запрограмми-
запрограммировано; откуда же может появиться случайность?
В этом вопросе и в самом деле есть некоторые труд-
трудности, но они относятся скорее к философии, так что мы
на них останавливаться не будем.
На всякий случай заметим только, что случайные
величины, о которых шла речь в § 2, — это идеальные
математические понятия. Вопрос о том, можно ли с их
помощью описать какое-либо явление природы, решает-
решается опытом. Такое описание всегда является приближен-
приближенным. Более того, случайная величина, которая вполне
удовлетворительно описывает какую-то физическую вели-
величину в одном круге явлений, может оказаться плохой
характеристикой этой же величины при исследовании
других явлений.
Точно так же дорога, которую на карте страны мож-
но считать прямой (идеальной математической прямой
«без ширины»), становится полосой с изгибами на круп-
крупномасштабном плане населенного пункта...
Обычно различают три способа получения случайных
величии: таблицы случайных чисел, генераторы случай-
случайных чисел и метод псевдослучайных чисел.
3.1. Таблицы случайных чисел. Проделаем
следующий опыт. Напишем на десяти одинаковых бу-
бумажках цифры 0, 1, 2, ..., 9. Положим эти бумажки в
шапку, перемешаем и будем извлекать оттуда по одной
бумажке, каждый раз возвращая ее назад и снова пере-
перемешивая все бумажки. Полученные таким образом циф-
цифры запишем в виде таблицы, подобной таблице А в конце
книги на стр. 62 (в таблице А цифры для удобства объ-
объединены в группы по 5 штук).
Такая таблица называется таблицей случайных чи-
чисел, хотя правильнее было бы назвать ее таблицей слу-
20

чайных цифр. Можно ввести ее в запоминающее устрой-
устройство ЭВМ. И в ходе расчета, когда нам понадобится
значение случайной величины с распределением
0 ' 2 ¦¦¦ ^	,22)
0,1 0,1 0,1
0,1
мы будем брать очередную цифру из этой таблицы.
Самая большая из опубликованных таблиц случай-
случайных чисел содержит 1 000 000 цифр*). Конечно, она со-
составлялась с помощью более современной техники, чем
шапка: была сконструирована специальная рулетка с
использованием электро-
электроники. Простейшая схема
такой рулетки приведена
на рис. 7 (вращающийся
диск резко останавливает-
останавливается и выбирается та циф-
цифра, на которую указывает
неподвижная стрелка).
Необходимо заметить,
что составить хорошую
таблицу случайных чисел
не так просто, как может
показаться. Любой реаль-
реальный физический прибор
вырабатывает случайные
величины с распределением, несколько отличающимся
от идеального распределения B2). К тому же в ходе
бпыта возможны ошибки (например, одна из бумажек на
некоторое время может пристать к подкладке). Поэтому
составленные таблицы тщательно проверяются с по-
помощью специальных статистических тестов: не противо-
противоречат ли те или иные свойства группы чисел гипотезе о
том, что эти числа — значения случайной величины B2).
Приведем один из простейших тестов. Рассмотрим
таблицу, содержащую N цифр. Пусть число нулей в этой
таблице v0, число единиц vi, число двоек v2 и т. д. Вы-
Вычислим сумму
Рис. 7.
*) RAND Corporation, A million random digits with 100 000 nor-
normal deviates, The Free Press, 1955.
21

Теория вероятностей позволяет предсказать, в каких
пределах может лежать эта сумма: она не должна быть
слишком большой (так как математическое ожидание
каждого из Vj равно 0,Ш), но не должна быть и слиш-
слишком малой (так как это означало бы «слишком законо-
закономерное» распределение значений).
Таблицы случайных чисел используются только при
расчетах по методу Монте-Карло вручную. Дело в том,
что все ЭВМ обладают сравнительно малой внутренней
памятью, и большая таблица туда не помещается. Хра-
Хранить же таблицу во внешней памяти и постоянно об-
обращаться к ней за числами — это очень замедляет
счет.
Однако не исключено, что со временем возможности
памяти у ЭВМ резко вырастут, и тогда таблицы слу-
случайных чисел станут весьма практичными.
3.2. Генераторы случайных чисел. Каза-
Казалось бы, что упоминавшуюся в 3.1 рулетку можно при-
присоединить к вычислительной машине и по мере надоб-
надобности вырабатывать случайные числа. Однако любой
механический прибор будет слишком медленным для
ЭВМ. Поэтому в качестве генераторов случайных вели-
величин чаще всего используют шумы в электронных лам*
пах: если за некоторый фиксированный промежуток вре-
времени At уровень шума превысил заданный порог четное
число раз, то записывается нуль, а если нечетное число
раз, то записывается единица*).
На первый взгляд это очень удобный способ. Пусть т
таких генераторов работают параллельно, работают все
время и засылают случайные нули и единицы во все
двоичные разряды специальной ячейки. Каждый такт —
одно m-разрядное число. В любой момент счета можно
обратиться к этой ячейке и взять оттуда значение слу-
случайной величины у. равномерно распределенной в интер-
интервале @, 1). Конечно, значение приближенное, записан-
записанное в форме т-разрядной двоичной дроби 0, a(i)ctB) •..
... Щт), где каждая из величин щ^ имитирует случайную
величину с распределением
О 1\
1/2 1/2/-
*) Существуют и более совершенные конструкции.
22

Однако н этот метод не свободен от недостатков. Во-
первых, трудно проверить «качество» вырабатываемых
чисел. Проверки приходится делать периодически, так
как из-за каких-либо неисправностей может возникнуть
так называемый «дрейф распределения» (то есть нули
и единицы в каком-либо из разрядов станут появляться
не одинаково часто). Во-вторых, обычно все расчеты на
ЭВМ проводят дважды, чтобы исключить возможность
случайного сбоя. Но воспроизвести те же случайные чис*
ла невозможно, если их по ходу счета не запоминать,
А если их запоминать, то мы снова приходим к сЛучаЮ
таблиц.
Датчики такого типа, несомненно, окажутся полез-*
ными тогда, когда будут строиться специализированные
ЭВМ для решения задач методом Монте-Карло. А для
универсальных ЭВМ, на которых расчеты с помощьк)
случайных чисел проводятся лишь изредка, содержать и
эксплуатировать специальное устройство просто не эко*
номично. Лучше использовать так называемые псевдо-
псевдослучайные числа.
3.3. Псевдослучайные числа. Поскольку «ка-
чество» используемых случайных чисел проверяется с
помощью специальных тестов, можно не интересоваться
тем, как эти числа получены: лишь бы они удовлетво-
удовлетворяли принятой системе тестов. Можно даже попытаться
вычислять их по заданной формуле. Но, конечно, это
должна быть весьма хитрая формула.
Числа, получаемые по какой-либо формуле и имити-
имитирующие значения случайной величины у, называются
псевдослучайными числами. Под словом «имитирующие»
подразумевается, что эти числа удовлетворяют ряду те-
тестов так, как если бы они были значениями этой случай-
случайной величины.
Первый алгоритм для получения псевдослучайных
чисел был предложен Дж. Нейманом. Он называется
методом середины квадратов. Поясним его на примере.
Пусть задано 4-значное целое число «i = 9876. Возве-
Возведем его в квадрат. Получим, вообще говоря, 8-значное
число /^ = 97 535 376. Выберем четыре средние цифры
из этого числа и обозначим я2=*5353.
Затем возведем п2 в квадрат (п\ = 28 654 609) и сно-
снова извлечем 4 средние цифры. Получим «3 = 6546.
23

Далее, я| = 42 850116, я4 = 8501; ftJ = 7
«5 = 2б70; «2 = 07128 900, я6=1289 и т. д.
В качестве значений у предлагалось использовать
значения 0,9876; 0,5353; 0,6546; 0,8501; 0,2670; 0,1289
и т. д.*).
Но этот алгоритм не оправдал себя: получалось
больше чем нужно малых значений. Поэтому разными
исследователями были разработаны другие алгоритмы.
Некоторые из них используют особенности конкретных
ЭВМ. В качестве примера рассмотрим один из таких
алгоритмов, применяемый на ЭВМ Стрела.
Пример**). Стрела — трехадресная ЭВМ с плавакдаей запя-
запятой. Ячейка, в которую записывается число х, состоит из 43 двоич-
двоичных разрядов (рис. 8). Машина работает с двоичными числами в
Т	Мантисса	I	Порядок
Знан мантисс/»	Зна/f порядно
Рис. 8.
нормализованной форме: x=±q-2*p, где р— порядок числа, q —<
мантисса. В /-м разряде ячейки может стоять нуль или единица;
обозначим эту величину е,. Тогда
Условие нормализации 0,5 <; <? < 1 означает, что 8i всегда равно 1.
Знак « + » изображается нулем, знак « — » —единицей.
Число v*.+i получается из числа ук тремя операциями:
1)	умножаем у* на большую константу, обычно 10";
2)	изображение произведения 1017 у* "сдвигаем на 7 разрядов
влево (так что первые 7 разрядов произведения пропадают, а в раз-
разрядах с 36-го по 42-й оказываются нули);
3)	берем абсолютную величину полученного числа (при этом
число нормализуется); это и будет "ул+ь
Если начинать с уо=1, то таким процессом получаются более
80 000 различных чисел у*, затем в последовательности возникает
*) Этот алгоритм можно записать в виде формулы
"M-i = f("*), где F означает совокупность операций, которые надо
проделать иад числом пк, чтобы получить пА+1. Число щ задается.
Псевдослучайные числа ул=10~4гсь
**) См. И. М. Соболь, Псевдослучайные числа для машины
Стрела, Теория верояти. и ее примен., 3, № 2, 1958, 205—21!.
24

период, и числа начинают повторяться. Различные проверки первых
50 000 чисел дали вполне удовлетворительные результаты. Эти числа
неоднократно использовались при решении самых разнообразных
задач.
Достоинства метода псевдослучайных чисел довольно
очевидны. Во-первых, на получение каждого числа за-
затрачивается всего несколько простых операций, так что
скорость генерирования случайных чисел имеет тот же
порядок, что скорость работы ЭВМ. Во-вторых, про-
программа занимает всего несколько ячеек накопителя,
В-третьих, любое из чисел yh может быть легко воспро*
изведено. В-четвертых, нужно лишь один раз проверить
«качество» такой последовательности, затем ее можно,
много раз безбоязненно использовать при расчете сход*
ных задач.
Единственный недостаток метода — ограниченность
«запаса» псевдослучайных чисел. Однако существуют
способы, позволяющие получать гораздо больше чисел,
В частности, можно менять начальные числа у0.
Подавляющее большинство расчетов по методу Мон-
Монте-Карло в настоящее время осуществляется с использо-
использованием псевдослучайных чисел.
§ 4. Преобразования случайных величин
При решении различных задач приходится моделиро-
моделировать различные случайные величины. На ранних этапах
использования метода Монте-Карло некоторые вычисли-
вычислители пытались для нахождения каждой случайной вели»
чины строить свою рулетку. Например, чтобы находить
значения случайной величины с распределением
0,5 0,25 0,125 0,125
можно использовать рулетку, изображенную на рис. 9,
действующую точно так же, как рулетка, изображенная
на рис. 7, но с неравномерными делениями, пропорцио-
пропорциональными р^
Однако это оказалось совершенно не нужным: зна-
значения любой случайной величины можно получить пу-
путем преобразования значений одной какой-либо (так
25

сказать «стандартной») случайной величины. Обычно
роль такой величины играет случайная величина у, равно-
равномерно распределенная в @, 1). Как получать значе-
значения у, мы уже знаем.
Условимся процесс нахождения значения какой-либо
случайной величины | путем преобразования одного или
нескольких значений ука-
указывать разыгрыванием
случайной величины |.
4.1. Р а з ыг р ы в а н и е
дискретной случай-
случайной величины. Допу-
Допустим, что нам нужно
получать значения слу-
случайной величины | с рас-
распределением
х2
Pi
Рп
Рис. 9.
Рассмотрим интервал
О< у < 1 и разобьем его
на п интервалов, длины которых равны ри рч, ..., рп.
Координатами точек деления, очевидно, будут у = р\,
У = Р\+Р2, У = Р\ + Р'2 + РЗ, -.., У=Р\+Р2+ ¦¦¦ +Рп-\-
Полученные интервалы занумеруем цифрами 1, 2,...
,.., п (рис. 10). На этом приготовления к розыгрышу |
Р,
Р/'Рг+Рз
1~Рп
Рис. 10.
заканчиваются. Каждый раз, когда нам надо будет «по-
«поставить опыт» и разыграть значение |, мы будем выби-
выбирать значение у и строить точку у = у. Если эта точка
попадет в интервал с номером I, то будем считать, что
\ = %i (в этом опыте).
Доказать законность такой процедуры совсем легко.
В самом деле, так как случайная величина у равномер-
равномерно распределена в @, 1), то вероятность того, что у ока-
26

жется в некотором интервале, равна длине этого интер»
вала. Значит,
P[0<V<pl]=Pl,
Р\Р\ <У<Р1 + Р2}=Р2>
Согласно нашей процедуре ? = *,- тогда, когда
а вероятность этого равна pt.
Конечно, на ЭВМ можно обойтись без рис. 10. Пред-
Предположим, что чиела*1, *2i •••, хп расположены в ячейках
¦
Находим
¦
V
Y < Р?
нет
Увеличиваем
на 1 адреса
ячеек, из кото-
которых берутся
Pi U X
i
да
Полагаем
\
Восстанавли-
Восстанавливаем адреса
Р\ и хх
Рас. 11.
накопителя подряд, а вероятности ри Pi+'P2, Pi + P2+:
,+ Рз. • • •> 1—тоже подряд. Блок-схема подпрограммы
для розыгрыша % приведена на рис. 11.
27

Пример. Разыграть 10 значений случайной величины 9 с рас-
распределением
6 ^ \0,58 0,42
В качестве значений у выберем пары цифр из таблицы А на
стр. 62, умноженные на 0,01 *). Итак, у=0,86; 0,51; 0,59; 0,07; 0,95;
0,66; 0,15; 0,56; 0,64; 0,34.
Очевидно, по нашей схеме значениям у меньшим, чем 0,58, от-
отвечают значения 6=3, а значениям у^0,ж — значения 6 = 4. Сле-
Следовательно, мы получим значения: 9 = 4; 3; 4; 3; 4; 4; 3; 3; 4; 3.
Заметим, что порядок нумерации значений Хи хг, ...
..., хп в распределении | может быть произвольным,
однако он должен быть фиксирован до начала розы-
розыгрыша.
4,2. Разыгрывание непрерывной случай-
случайной величины. Допустим теперь, что нам нужно по-
получить значения случайной величины %, распределенной
В интервале (а, Ь) с плотностью р(х).
Докажем, что значения % можно находить из урав-
уравнения
I
p(x)dx=y,	B3)
то есть выбрав очередное значение у, надо решить урав-
уравнение B3) и найти очередное значение |.
Для двказательства рассмотрим функцию (рис. 12)
х
у = j P(x)dx.
а
Из общих свойств плотности A5) и A6) следует, что
у(а) = 0, у(Ь)=1,
а производная
*) Так как в этом примере р, заданы с двумя десятичными
знаками, то вполне достаточно брать значения y с двумя десятич-
десятичными знаками. При таком приближении возможен случай у=0,58,
который надо объединить со случаем y>0,58 (ибо значение у=0,00
возможно, а значение у='.00 невозможно). При использовании
многозначных у случай равенства y — pi мало вероятен и его можно
относить к любому из неравенств.
28

Значит, функция у(х) монотонно возрастает от 0до 1.
И любая прямая у —у, где 0< у < 1, пересекает график
у=у(х) в одной-единственной точке, абсциссу которой
мы и принимаем за |. Таким образом, уравнение B3)
всегда имеет одно и только одно решение.
Выберем теперь произвольный интервал (а', Ь'), со-
содержащийся внутри (а,Ь). Точкам этого интервала
а' < х < Ь'
отвечают ординаты кривой у = у(х), удовлетворяющие
неравенству
у{а')<У<у(Ь').
Поэтому если | принадлежит интервалу а' < х < Ь', то
у принадлежит интервалу у(а')<у <y(bf), и наоборот
(рис. 13). Значит,
Р {а' < |< Ь'}=Р{у(а') < у <у{Ь')}.
Так как y равномерно распределена в @, 1), то
Ь'
Р [у (а1) <у<у (&')} == У\Ь') - U Ю =
а'
Итак,
J p{x)dx.
Р{а' <\<Ь'] = \ p(x)dx,
а это как раз и означает, что случайная величина ?, яв-»
ляющаяся корнем уравнения B3), имеет плотность
вероятностей р(х).
29

Пример. Случайная величина г| называется равномерно рас-
пределенной в интервале (а,Ь), если ее плотность постоянна в этом
интервале:
р {х) = -т	 при а < х < Ь.
.Чтобы разыгрывать значения г|, составим уравнение B3):
11
I
dx
Интеграл легко вычисляется:
Ъ~а
= Y-
Отсюда получается явное выражение для г|:
ц = а-\-ч(Ь — а).	B4)
Другие примеры применения формулы B3) приведе-
приведены в пп. 5.2 и 8.3.
4.3. Метод Неймана для разыгрывания не-
непрерывной случайной величины. Может ока-
оказаться, что разрешить уравнение B3) относительно g
весьма трудно. Напри-
Например, в случае, когда инте-
интеграл от р(х) не выра-
выражается через элементар-
элементарные функции или когда
плотность р(х) задана
графически.
Предположим, что слу-
случайная величина | опре-
определена на конечном интер-
интервале (а, Ь) и плотность ее
ограничена (рис. 14):
Рис- 14-
Разыгрывать значение |
можно следующим обра-
образом:
чины
тами
зом:
1) выбираем два значения у' и у" случайной вели-
велиы y и строим случайную точку Г (г/; rj") с координа-
координаи '
30

2) если точка Г лежит под кривой у — р(х), то пола-
полагаем ?, = ц', если же точка Г лежит над кривой у = р(х),
то пару (у', у") отбрасываем и выбираем новую пару
значений (у', у").
Обоснование этого метода приведено в п. 9.1.
4.4. О разыгрывании нормальных вели-
величин. Существует много весьма разнообразных приемов
для разыгрывания различных случайных величин. Мы
не будем здесь на них останавливаться. Они исполь-
используются обычно лишь тогда, когда приемы пп. 4.2 и 4.3
оказываются малоэффективными.
В частности, такой случай имеет место для нормаль-
нормальной случайной величины ?, так как уравнение
в явном виде неразрешимо, а интервал возможных зна-
значений ? бесконечен.
В таблице Б в конце книги (на стр. 62) приведены
значения (уже разыгранные) нормальной случайной ве-
величины % с математическим ожиданием М? = 0 и с дис-
дисперсией D?=l. Нетрудно доказать*), что случайная
величина
?' = а + о?	B5)
будет также нормальной, причем из A0) и A1) сле-
следует, что
D?' 2
Таким образом, формула B5) позволяет с помощью
таблицы Б разыгрывать значения любых нормальных
величин.
4.5. Снова о примере из п. 1.2. Теперь можно
объяснить, как выбирались случайные точки на рнс. 1 и 2.
На рис. 1 построены точки с координатами
Значения у' и у" вычислялись по пятеркам цифр из таб-
таблицы А: *! = 0,86515; г/, = 0,90795; л-2 = 0,66155; г/2 = 0,66434
и т. д.
*) Доказательство приведено в п. 9.2.
31

Можно доказать*), что так как абсциссы и орди-
ординаты наших точек независимы, то вероятность попада-
попадания такой точки в любую область внутри квадрата рав-
равна площади этой области. Иначе говоря, точки равно-
равномерно распределены в квадрате.
На рис. 2 построены точки с координатами
причем значения tf и Z," выбирались из таблицы Б под-
подряд:
Х! = 0,5 + 0,2-0,2005, #1 = 0,5 + 0,2- 1,1922;
х, = 0,5 + 0,2 (—0,0077),...
Одна из точек, оказавшаяся вне квадрата, отброшена.
Из формулы B5) следует, что абсциссы и ординаты
этих точек представляют собой нормальные случайные
величины со средними о = 0,5 и с дисперсиями а2 = 0,04.
*) Доказательство приведено в п. 9.3.

ГЛАВА 2
ПРИМЕРЫ ПРИМЕНЕНИЯ МЕТОДА МОНТЕ-КАРЛО
§ 5. Расчет системы массового обслуживания
5.1.	Описание задачи. Рассмотрим одну из са-
самых простых систем массового обслуживания. Система
эта состоит из п линий (или каналов, или пунктов об-
обслуживания), каждый из которых может «обслуживать
потребителей». В систему поступают заявки, причем мо-
моменты их поступления случайные. Каждая заявка посту-
поступает на линию номер 1. Если в момент поступления k-w.
заявки (назовем его Th) эта линия свободна, то она при-
приступает к обслуживанию заявки, что продолжается t3 ми-
минут (t3— время занятости линии). Если в момент 7\ ли-
линия номер 1 занята, то заявка мгновенно передается на
линию номер 2. И так далее...
Наконец, если все п линий в момент Th заняты, то си-
система выдает отказ.
Требуется определить, сколько (в среднем) заявок
обслужит система за время Т и сколько отказов она
даст?
Ясно, что задачи такого типа встречаются при иссле-
исследовании организации работы любых предприятий, а не
только предприятий бытового обслуживания. В некото-
некоторых очень частных случаях удается найти аналитиче-
аналитические решения. Однако в сложных случаях (о них будет
сказано ниже) метод Монте-Карло оказывается един-
единственным методом расчета.
5.2.	Простейший поток заявок. Первый во-
вопрос, возникающий при рассмотрении такой системы:
что представляет собой поток поступающих заявок? Этот
вопрос решается опытом, путем достаточно длительного
3 И. М. Соболь	33

наблюдения за заявками. Изучение потоков заявок в раз-
различных условиях позволило выделить некоторые доста-
достаточно часто встречающиеся случаи.
Простейшим потоком (или потоком Пуассона) назы-
называется такой поток заявок, когда промежуток времени х
между двумя последователь-
последовательными заявками есть случай-
случайная величина, распределен-
распределенная в интервале @, оо) с
плотностью
р(х) — ае~ах. B6)
Закон B6) называют также
экспоненциальным распре-
х делением (см. рис. 15, где
построены плотности B6)
при а= 1 и о = 2).
Легко вычислить математическое ожидание х:
оо	оо
Мт = J xp(x)dx— f xae~axdx.
После интегрирования по частям {и-х, dv = ae~axdx)
получим
Мт = [—.
1
Параметр а называется плотностью потока заявок.
Формулу для розыгрыша т легко получить из уравне-
уравнения B3), которое в нашем случае запишется так:
Вычислив интеграл, стоящий слева, получим соотношение
откуда
34

. Впрочем, величина 1—у распределена точно так же,
как y> и поэтому можно вместо последней формулы
использовать формулу
T==_ilnY.	B7)
5.3. Схема расчета. Итак, рассмотрим работу си-*
стемы из п. 5.1 в случае простейшего потока заявок.
Каждой линии поставим в соответствие ячейку вну-
внутреннего накопителя ЭВМ, в которую будем записывать
момент, когда эта линия освобождается. Обозначим мо-
момент освобождения i-й линии через ti. За начальный мо*
мент расчета выберем момент поступления первой заяв-
заявки 7 = О. Можно считать, что в этот момент все tt равны
Т\: все линии свободны. Время окончания расчета
Ткон— Т\ + Т.
Первая заявка поступает на линию номер 1. Значит,
в течение t3 линия эта будет занята. Поэтому мы должны
заменить t\ на новое значение (^i)Hob=7'i + 4, добавить
единицу к счетчику выполненных заявок и перейти к
рассмотрению второй заявки.
Предположим, что k заявок уже рассмотрены. Тогда
надо разыграть момент поступления (k + \)-\\ заявки.
Для этого выбираем очередное значение у и по фор-
формуле B7) вычисляем очередное значение т = т^. А затем
вычисляем момент поступления
Свободна ли в этот момент первая линия? Для уста-
установления этого надо проверить условие
'i<W	B8)
Если это условие выполнено, то, значит, к моменту 7",i+t
линия уже освободилась и может обслужить эту заявку.
Мы должны заменить t\ на Tk+i + t3, добавить единицу
к счетчику выполненных заявок и перейти к следую-
следующей заявке.
Если условие B8) не выполнено, то это значит, что
первая линия в момент Th+i занята. Тогда мы проверяем,
свободна ли вторая линия:
B<Tk+l?	B9)
3*	35

Если условие B9) выполнено, то мы заменяем t2 на
ts, добавляем единицу к счетчику выполненных
заявок и переходим к следующей заявке.
Если же условие B9) не выполнено, то переходим к
проверке условия
Может оказаться, что при всех / от 1 до п
то есть все линии в момент Tk+l заняты. Тогда надо до-
добавить единицу в счетчик отказов и потом перейти к рас-
рассмотрению следующей заявки.
Каждый раз, вычислив Тк+\, надо проверять еще ус-
условие окончания опыта
' k+l > ' кон-
Когда это условие окажется выполненным, опыт закан->
чивается. В счетчике выполненных заявок и в счетчике
отказов будут стоять числа цвып и цОтк.
Такой опыт повторяется N раз (с использованием
различных у)- И результаты всех опытов осредняются:
N
Л
^1°тк О')'
где (Лвыпу) и Hotk(j) — значения цвып и [i0Ti;, полученные в
/-м опыте.
На рис. 16 приведена блок-схема программы, осуще-
осуществляющей такой расчет (в случае необходимости мож-
можно получить значения цВыпо) н Мотко) для отдельных
испытаний в блоке «конец опыта»).
5.4. Более сложные задачи. Очень легко пока*
зать, что этот же метод позволяет рассчитывать несрав-
несравненно более сложные системы. Например, величина /3
может быть не постоянной, а случайной и различной
для различных линий (что соответствует различному
оборудованию или различной квалификации обслужи-
обслуживающего персонала). Схема расчета останется в основном
36

Ввод
j
данных
i
= 1
Конец опыта
нет
/нов — Уст Т
[да
da
да
Добавляем 1
к счетчику
выполненных
заявок
Подсчет резуль-
результатов, выдача,
конец задачи
да
.нет
'нов—'ст~Ь'
Добавляем I
к счетчику
отказов
1
т* = -
i
1
h 	 h
"¦нов — K
Iln
a
» +
\
Y
т*
1
Рис. 16.
37

такой же, но значения t3 придется каждый раз разыгры-
разыгрывать, и формула разыгрывания для каждой линии будет
своя.
Можно рассматривать так называемые системы с
ожиданием, в которых отказ выдается не сразу: заявка
хранится в течение некоторого времени ta (время пре-
пребывания заявки в системе), и если за это время какая-
нибудь линия освободится, то она обслужит эту заявку.
Можно рассматривать системы, в которых очеред-
очередную заявку принимает та линия, которая раньше всех
освободилась. Можно учесть случайный выход из строя
отдельных линий и случайное время ремонта каждой из
них. Можно учесть изменения плотности потока заявок
во времени. И многое, многое другое.
Конечно, расчеты таких систем не даются даром.
Чтобы получить результаты, имеющие практическую цен-
ценность, надо выбрать хорошую модель. Для этого прихо-
приходится достаточно тщательно изучать действительные по-
потоки заявок, проводить хронометраж работы отдельных
узлов и т. п.
Вообще надо знать вероятностные законы функцио-
функционирования отдельных частей системы. Тогда метод
Монте-Карло позволяет вычислить вероятностные зако-
законы работы всей системы, как бы сложна она ни была.
Такие методы расчета чрезвычайно полезны при
планировании предприятий: вместо дорогостоящего (а
иногда просто невозможного) эксперимента в натуре мы
можем экспериментировать на ЭВМ, моделируя разные
варианты организации работы пли использования обо-
оборудования.
§ 6. Расчет качества и надежности изделий
6.1. Простейшая схема расчета качеств а.
Рассмотрим изделие S, состоящее из некоторого (может
быть, большого) числа элементов. Например, если S
представляет собой электрический прибор, то элемента-
элементами его могут быть сопротивления {R(h)), емкости (C(ft)),
лампы и т. п. Предположим, что качество изделия опре-
определяется значением одного выходного параметра U, ко-
которое можно вычислить, зная параметры всех элементов
, Ябу, ..-; СA), СB), ...; ...). C0)
38

Если, например, U — это напряжение на рабочем
участке электрический цепи, то но законам Ома можно
составить уравнения для цепи и, решая их, найти U.
В действительности, однако, параметры элементов не
равны в точности указанным значениям. Например, со-
сопротивление, изображенное на рис. 17, может оказаться
любым в интервале от 20,9 до
23,1 /ш.	Гг—щр	П
Возникает вопрос: как повлияют [_] 22KS25%2W N
отклонения параметров всех эле-
элементов от номинальных на значе-	Рис. 17.
ние U?
Можно попытаться оценить пределы изменения U4
выбирая для всех элементов «худшие» значения пара-
параметров. Однако далеко не всегда известно, какой набор
параметров будет «худшим». К тому же, если число эле-
элементов велико, то такая оценка окажется сильно завы-
завышенной: на самом деле маловероятно, чтобы все пара-
параметры одновременно оказались наихудшими.
Поэтому разумнее считать параметры всех элементов
и саму величину U случайными величинами и попытать-
попытаться оценить математическое ожидание MU и дисперсию
DU. Величина MU — это среднее значение U для всей
партии изделий, а величина DU показывает, какие от-
отклонения U от MU будут встречаться на практике.
Напомним (об этом было сказано в н. 2.2), что
MU^f(MR{l), MR{2), ...; МСA), МС,2), ...; ...).
Вычислить аналитически распределение U при мало-
мальски сложной функции f невозможно. Иногда это
можно сделать экспериментально, просмотрев большую
партию готовых изделий. Но это возможно далеко не
всегда, а на стадии проектирования — никогда.
Попробуем применить метод Монте-Карло. Для этого
нужно: а) знать вероятностные характеристики всех эле-
элементов, б) знать функцию f (точнее, уметь вычислять
значение U по любым фиксированным значениям /?(i>,
R{2), ¦ ¦ -\ C(i), СB), ...; ...).
Вероятностное распределение параметров для каждо-
каждого отдельного элемента можно получить эксперимен-
экспериментально, путем просмотра большой партии таких эле-
элементов. Весьма часто это распределение оказывается
нормальным. Поэтому многие исследователи поступают
39

следующим образом: например, считают сопротивление
элемента, изображенного на рис. 17, нормальной слу-
случайной величиной р с математическим ожиданием Мр=
=22 и с За== 1,1 (напомним, что получить в одном опыте
значение р, отклоняющееся от Мр больше чем на За,
практически невозможно, см. B0)).
Схема расчета оказывается очень простой: для ка«
ждого элемента разыгрывается значение параметра; за-
затем по формуле C0) вычисляется значение U. Повторив
этот опыт N раз и получив значения Uu U2, ..., UN, мо-
можем приближенно считать, что
j-l
DU
_!_Гу«,-_1(уИ1
/
—
—
—
4
При больших N в последней формуле можно заме*
нить множитель l/(N—1) на 1/N, и тогда эта формула
окажется простым следствием формул (8) и (9). В ма-
математической статистике доказывается, что при неболь-
небольших N лучше сохранить множитель l/(N—1).
6.2. Примеры расчета надежности. Пусть
мы хотим оценить среднее время безотказной работы
изделия, предполагая, что
известны характеристики
безотказной работы ка-
каждого из элементов.
Если считать, что вре-
Рис. 18.	мя безотказной работы
каждого элемента t@)—¦
фиксированная величина, то вычислить время безотказ-
безотказной работы t изделия не составит труда. Например, для
изделия, схематически изображенного на рис. 18, в ко-
котором выход из строя одного элемента влечет за собой
выход из строя всего изделия,
А для изделия, схематически изображенного на рис. 19,
в котором один из элементов дублирован,
* = 1Шп[/A); ti2y, тах(^з); tw); ti5)\,	C2)
40

так как если, например, элемент № 3 выйдет из строя,
то изделие будет продолжать работать на одном эле-
элементе № 4.
В действительности время безотказной работы лю-
любого элемента представляет собой случайную величину
в(й). Когда мы говорим, что срок службы электрической
лампочки равен 1000 часов, то это лишь среднее значе-
значение Мв величины в: всем
известно, что одна лам-
лампочка перегорает быстрее,
а другая (в точности та-
такая же) горит дольше.
Если известны плотно-
плотности распределений в^)для
каждого из элементов из-
изделия, то Мв можно со-
сосчитать методом Монте-Карло в точности так же, как
это делалось в п. 6.1. В самом деле, для каждого эле-
элемента можно разыграть значение величины в^), пусть
это будет t{k)- Затем по соответствующей формуле C1)
или C2) можно вычислить значение t случайной вели-
величины в. Повторив этот опыт достаточно много (N) раз,
можем считать, что
N
./
/
—
2
3
4
4
\
/
j
Рис. 19.
где tj — значение t, полученное в /-м опыте.
Необходимо заметить, что вопрос о распределении
времени жизни ©^ для отдельных элементов не так уж
прост: для наиболее долговечных элементов организа-
цня эксперимента затруднительна, так как нужно до-
дождаться, пока выйдет из строя достаточно много эле-
элементов.
6.3. Дальнейшие возможности метода.
Приведенные примеры показывают, что методика расче*
та качества проектируемых изделий весьма проста по
идее. Нужно знать вероятностные характеристики всех
элементов изделия и уметь вычислить интересующую нас
величину, как функцию от параметров этих элементов.
Тогда случайность параметров можно учесть путем мо-
моделирования.
41

При моделировании можно получить гораздо боль-
больше полезной информации, а не только математическое
ожидание и дисперсию интересующей нас величины.
Пусть, например, мы получили большое число N значе-
значений U\, U2, ¦ ¦., UN случайной величины U. По этим зна-
значениям можно построить приближенную плотность рас-
распределения U. Вопрос этот относится, по существу, к
статистике, так как речь идет об обработке результатов
експериментов (только проведены они на ЭВМ). Поэтому
ограничимся лишь конкретным примером.
Допустим, что мы получили всего iV= 120 значений
Uи U2, ¦ ¦ -, Una случайной величины U, и все они заклю-
заключены в пределах
l<Uj< 6,5.
Разобьем интервал 1 < х < 6,5 на 11 (любое число,
не слишком большое и не слишком малое) равных ин-
интервалов длиной Дх = 0,5 и сосчитаем, сколько значе-
значений Uj попали в каждый интервал. Числа попаданий
приведены на рис. 20.
2 О I 15 24 32 /7 19 7 t 2
_i	1	1	1	1	1	(	1	1	1	1	1	p_».
/	2	3	4	5	6	7 X
Рис. 20.
Частота попадания в какой-либо интервал получается
делением числа попаданий на N= 120. В нашем примере
частоты равны: 0,017; 0; 0,008; 0,12; 0,20; 0,27; 0,14; 0,16;
0,06; 0,008; 0,017.
Над каждым из интервалов разбиения построим
прямоугольник, площадь которого равна частоте попа-
попадания Uj в этот интервал (рис. 21). Иначе говоря, вы-
высота каждого прямоугольника равна частоте, деленной
на Аде. Полученную ступенчатую линию называют гисто-
гистограммой.
Гистограмма служит приближением к неизвестной
плотности случайной величины U. Поэтому, например,
площадь гистограммы, заключенная между х = 2,5 и
х = 5,5, дает нам приближенное значение вероятности
Р {2,5 <U< 5,5} ^0,94.
Следовательно, на основании проведенного расчета
(опыта) можно считать, что с вероятностью, приблнзн-
42

тельно равной 0,94, величина U заключена в интервале
2,5 <U< 5,5.
На рис. 21 для сравнения построена плотность нор«
мальной случайной величины tf с параметрами а = 3,85,
о=0,88. Если по этой плотности вычислить вероятность
того, что ?' заключена в интервале 2,5 < ?' < 5,5, то по«
лучим довольно близкое значение 0,91 *).
*) Ход вычислений. Согласно A4) запишем
5,5
(х-аУ
dx.
В интеграле сделаем замену переменной x—a = at. Тогда получим
Р {2,5
5,5) =
1
2к
dt.
, 2,5 — а	, ,. , 5,5 — a QQ
где ^ =	= — 1,54, г2 =	= 1,оо. Значение послед-
последнего интеграла вычисляется с помощью таблиц так называемого ин-
интеграла вероятностей Ф(х), в которых приведены значения функции
х р
Ф (л-) = ^Lr [ e" dt.
Получим
Р {2,5 < $' < 5,5} = 0,5 [Ф A,88) + Ф A,54)] = 0,91.
43

6.4. Замечание. К сожалению, расчеты такого ти-
типа пока проводят редко. В чем главная причина этого —
трудно сказать. Скорее всего в том, что конструкторы н
проектировщики не знают о такой возможности.
К тому же, прежде чем рассчитывать таким методом
какие-либо изделия, надо изучить вероятностные харак-
характеристики всех элементов, входящих в эти изделия. Это
немалая работа. Правда, зная эти характеристики, мож-
можно оценивать качество любых изделий, состоящих из
этих элементов. Можно оценивать изменение качества
при замене одних элементов другими.
Нужно надеяться, что в ближайшие годы такие рас-
расчеты станут привычным делом. А вероятностные харак-
характеристики элементов будут всегда выдаваться предприя-
предприятиями, их изготовляющими.
§ 7. Расчет прохождения нейтронов сквозь
пластинку
Вероятностные законы взаимодействия отдельной эле-
элементарной частицы (нейтрона, фотона, мезона и др.)
с веществом известны. Обычно требуется найти макро-
макроскопические характеристики процессов, в которых уча-
участвует огромное количество таких частиц: плотности,
потоки и т. п. Эта ситуация похожа на ту, с которой мы
сталкивались в §§ 5 и 6, и весьма удобна для использо-
использования метода Монте-Карло.
Пожалуй, чаще всего метод Монте-Карло исполь-
используется в нейтронной физике. Мы рассмотрим простей-
простейший вариант задачи о прохождении нейтронов сквозь
пластинку.
7.1. Постановка задачи. Пусть на однородную
бесконечную пластинку 0^.x4^.h падает поток нейтро-
нейтронов с энергией Ео. Угол падения 90°. При столкновении
с атомами вещества, из которого состоит пластинка,
нейтроны могут упруго рассеиваться или поглощаться.
Предположим для простоты, что энергия нейтрона при
рассеянии не меняется и любое направление «отскока»
нейтрона от атома одинаково вероятно (последнее иногда
справедливо в веществах с тяжелыми атомами). На
рис. 22 изображены различные судьбы нейтронов: нейт-
нейтрон (а) проходит сквозь пластинку, нейтрон (б) погло-
поглощается, нейтрон (в) отражается пластинкой.
44

Требуется вычислить вероятность прохождения нейт-
нейтрона сквозь пластинку р+, вероятность отражения ней-
нейтрона пластинкой р~ и вероятность поглощения нейтрона
в пластинке р°.
Взаимодействие нейтронов с веществом характери-
характеризуется в рассматриваемом случае двумя постоянными
2С и Е8( которые называют-
называются сечением поглощения и
сечением рассеяния. Индек-
Индексы с и s — это первые буквы
английских слов capture
(захват) и scattering (рас-
(рассеяние).
Сумма этих сечений на-
называется полным сечением
2 = 2,
Рис. 22.
Физический смысл сече-
сечений следующий: при стол-
столкновении нейтрона с атомом
вещества вероятность погло-
поглощения равна 2с/2, а вероят-
вероятность рассеяния равна 2.,/Е.
Длина свободного пробега нейтрона К (то есть длина
пути от столкновения до столкновения) — это случайная
величина. Она может принимать любые положительные
значения с плотностью вероятностей
Нетрудно заметить, что эта плотность величины X со-
совпадает с плотностью B6) случайной величины т для
простейшего потока заявок. По аналогии с п. 5.2 мы мо-
можем сразу написать выражение для средней длины сво-
свободного пробега
МЯ, = 1/2
и формулу для разыгрывания X:
Остается еще выяснить, как выбирать случайное на-
направление нейтрона после рассеяния. Так как задача
симметрична относительно оси х, то направление вполне
45

определяется одним углом ф между направлением ско-
скорости нейтрона и осью Ох. Л1ожно доказать*), что тре-
требование равной вероятности любого направления в этом
случае равносильно требованию, чтобы косинус этого
угла ц = соБф был равномерно распределен в интервале
(—1, 1). Из формулы B4) при а = — 1, Ь=\ следует фор-
формула для разыгрывания ц:
7.2. Схема расчета путем моделирования
истинных траекторий. Предположим, что нейтрон
испытал k-e рассеяние
h т.
внутри пластинки в точке с
абсциссой xh и после этого на-
начал двигаться в направле-
направлении Lift.
Разыграем длину свободно-
свободного пробега
Л, = -A/2) In Y
и вычислим абсциссу следую-
следующего столкновения (рис. 23)
Рис. 23.
xk + l —-X^-
Проверим условие прохо-
прохождения сквозь пластинку:
Если это условие выполнено, то счет траектории нейтро-
нейтрона заканчивается и добавляется единица к счетчику про-
прошедших частиц. В противном случае проверяем условие
отражения:
xk+1 < 0.
Если это условие выполнено, то счет траектории закан-
заканчивается и добавляется единица к счетчику отраженных
частиц. Если же и это условие не выполнено, то есть
0 ^xh+\^Ch, то, значит, нейтрон испытал (?+1)-е стол-
столкновение внутри пластинки, и надо разыграть «судьбу»
нейтрона при столкновении.
*) Доказательство приведено в п. 9.4.
46

Согласно п. 4.1 выбираем очередное значение y и про*
веряем условие поглощения:
Y < 2с/2.
Если последнее неравенство выполнено, то счет траекто-
траектории заканчивается и добавляется единица к счетчику по-
поглощенных частиц. В противном случае мы считаем, что
нейтрон испытал рассеяние в точке с абсциссой xh+i.
Тогда разыгрываем новое направление скорости ней-
нейтрона
и затем повторяем весь цикл снова (но, конечно, уже с
другими значениями у).
Все y написаны без индексов, так как имеется в виду,
что каждое значение y используется всего один раз. Для
расчета одного звена траектории нужны три значения y-
Начальные значения для каждой траектории:
После того как будут сосчитаны N траекторий, ока-
окажется, что N+ нейтронов прошли сквозь пластинку, N~
нейтронов отразились от нее, а № нейтронов были по-
поглощены. Очевидно, искомые вероятности приближенно
равны отношениям
N+	ЛГ	„ Л'0
Р ~ N • Р ~~ N ' Р ^ N •
На рис. 24 приведена блок-схема программы для рас-
расчета этой задачи. Индекс / — это номер траектории, ин-
индекс k — номер столкновения (вдоль траектории).
Такая методика расчета хотя и очень естественна, но
несовершенна. В частности, таким методом трудно вы-
вычислить вероятность р+, когда она очень мала. А с таким
случаем как раз приходится сталкиваться при расчете
защиты от излучений.
Существуют более «хитрые» разновидности метода
Монте-Карло, позволяющие рассчитывать и такие слу-
случаи. Остановимся коротко на одном из простейших ва-
вариантов расчета с помощью так называемых «весов».
7.3. Схема расчета с использованием ве-
весов, заменяющих поглощение. Рассмотрим ту
же задачу о прохождении нейтронов. Предположим, что
47

HOB
л'н-ов=
№ -
"нов"
]f у
Jaoi
Уст
Juoi
да
Ввод данных
У
1
СТ
у
+1
у
нет
1
'¦ft у
•V,+ ,=xH
У
У
- In Y
Л?
нет
**+! <0?
da
da
ВыадсЛ0яиг
нет
г
¦ 1
V-H +
9
1 — -V
— 1
У
"НОВ =г= «СТ ~f~ ^
I
, Р°
Выдача результатов, конец задачи
Рис. 24.
48

вдоль одной и той же траектории движется «пакет», со-
состоящий из большого числа и>0 одинаковых нейтронов.
При столкновении в точке с абсциссой Xi количество по-
поглощенных нейтронов из «пакета» в среднем равно
шоBс/2), а количество нейтронов, испытавших рассея-
рассеяние, в среднем равно шоB8/2).
Добавим величину дооBс/2) к счетчику поглощенных
частиц, а за движением рассеянного «пакета» будем сле-
следить дальше, предположив, что весь оставшийся «пакет»
рассеялся в одном направлении.
Все формулы счета, приведенные в п. 7.2, остаются
прежними. Только при каждом столкновении количество
нейтронов в «пакете» будет уменьшаться
так как часть «пакета», содержащая Wk(Sci) нейтро-
нейтронов, будет поглощаться. И траектория теперь не может
закончиться поглощением.
Величину Wh обычно называют весом нейтрона и
вместо того, чтобы говорить о «пакете», состоящем нз
Wu нейтронов, говорят об одном нейтроне с весом Wh-
Начальный вес а>о обычно полагают равным 1. Это не
противоречит нашим рассуждениям о «большом пакете»,
ибо нетрудно заметить, что все до*, получающиеся при
расчете одной траектории, содержат w0 общим множн-
телем.
Блок-схема программы для реализации такого рас-
расчета приведена на рие. 25. Очевидно, она ничуть не слож-
сложнее, чем блок-схема, изображенная на рнс. 24. Однако
можно доказать*), что расчет р+ методом настоящего
пункта всегда выгоднее, чем расчет по методу п. 7.2.
7.4. Замечание.Существует довольно много разно-
разнообразных способов расчета, использующих различные
веса. Мы не можем здесь на них останавливаться.
Отметим только, что метод Монте-Карло позволяет
решать гораздо более сложные задачи об элементарных
частицах: исследуемая среда может состоять из различ-
различных веществ и иметь любую геометрическую структуру;
энергия частицы при каждом столкновении может
меняться. Можно учитывать много других ядерных
*) Доказательство приведено в п. 9.5.
4 и. м. Соболь	49

Ввод данных
У==1
0, ха = О, ц0 == 1, w0 = 1
4-
Ф
Ф
Ф
- In Y
VXkV-k
A?
HOB
*;г»=*с~т-н»*
da
нет
4-
л> + ,<0?
da
<
нет
Унов — Уст
ф
+ 1
УИов<ЛГ?
da
нет
	V
= 2v—1
Ф
^НОВ ~~ '
Вычисление р+, р~,
Ф
Выдача результатов,
конец задачи
Рис. 25.
50

процессов (например, возможность деления атома при
столкновении с нейтроном и образование при этом но-
новых нейтронов). Можно рассчитать условие возникнове-
возникновения и поддерживания цепной реакции и т. д. и т. п.
§ 8. Вычисление определенного интеграла
Задачи, рассмотренные в §§ 5, 6, 7, были вероятност-
вероятностными по своей природе, и использование метода Монте-
Карло для их решения казалось довольно естественным.
Здесь рассматривается просто математическая задача:
приближенное вычисление определенного интеграла.
Так как вычисление определенных интегралов равно-
равносильно вычислению площадей, то можно было бы ис-
использовать метод примера 1.2. Вместо этого мы изложим
другой, более эффективный метод, позволяющий строить
разные вероятностные модели для решения этой задачи
методом Монте-Карло, и укажем, как среди них выби-
выбирать модели «получше».
8.1. Метод вычисления. Рассмотрим функцию
g(x), заданную на интервале а < х < Ь. Требуется при-
приближенно вычислить интеграл
ь
I=jg(x)dx.	C3)
Выберем произвольную плотность распределения
р\{х), определенную на интервале (а, Ь) (то есть произ-
произвольную функцию р\{х), удовлетворяющую условиям
A5) и A6)).
Наряду со случайной величиной |, определенной в ин-
интервале (о, Ь) с плотностью ръ(х), нам понадобится слу-
случайная величина
1\ = е AI Ръ (&)•
Согласно A8)
ь
Мл =j[g {х)!ръ (х)] р1 (х) dx = /.
Рассмотрим теперь iV одинаковых случайных величин
¦Ль тJ> •¦¦. 1» и применим к их сумме центральную
4*	51

предельную теорему п. 2.4. Формула B1) в этом случае
запишется так:
f M
1 Ж'..	, . о,/ «Ч !	0997	^34)
Последнее соотношение означает, что если мы выбе-
выберем N значений |ь |2, ..., lN, то при достаточно боль-
большом N
N
C5)
Оно показывает также, что с очень большой вероятно-
вероятностью ошибка приближения C5) не превосходит 3 YO^/N*
8.2. Как выбирать схему расчета. Мы ви-
видели, что для расчета интеграла C3) можно использо-
использовать любую случайную величину |, определенную в ин-
интервале (а, Ь). В любом случае
Однако дисперсия Dr], а с ней и оценка погрешности
формулы C5) зависят от того, какую величину % мы
используем. В самом деле,
Dti = М (л2) - Р = J [g2 (x)/Pi (х)] dx - Р.
а
Можно доказать*), что это выражение будет мини-
минимальным тогда, когда р\(х) пропорциональна \g{x)\.
Конечно, выбирать очень сложные р\{х) нельзя, так
как процедура разыгрывания значений | станет очень
трудоемкой. Но руководствоваться этой рекомендацией
при выборе р\(х) можно (см. пример п. 8.3).
На практике интегралы вида C3) методом Монте-
Карло не вычисляют: для этого есть более точные мето-
методы — квадратурные формулы. Однако при переходе к
многократным интегралам положение меняется: квадра-
квадратурные формулы становятся очень сложными, а метод
Монте-Карло остается почти без изменений.
*) Доказательство приведено в п. 9.6.
52

8.3. Ч и с л е н н ы й пример. Вычислим приближенно интеграл
я/2
/ = sin х dx.
о
Точное значение этого интеграла известно:
Я/2
f sin х dx = [— cos x]q/2 = 1.
О
Мы используем для вычисления две различные случайные вели-
величины %: с постоянной плотностью pi(x)~ 2/я (то есть | равномерно
распределена в интервале @, я/2)) и с линейной плотностью
Рис. 26.
Pl(x) =8*/я2. Обе эти плотности вместе с подынтегральной функ-
функцией sin х построены на рис. 26. Из этого рисунка видно, что линей-
линейная плотность более соответствует рекомендации п. 8.2 о том, что
желательна пропорциональность р| (х) и sin x. Поэтому нужно ожи-
ожидать, что второй способ вычисления даст лучший результат.
(а) Пусть р|(*)=2/я на интервале @, я/2)." Формула для
() Пусть р|(*J/я на интервал
разыгрывания \ может быть получена из
и -6 = я/2:
(, /)	ру
формулы B4) при
t =
-Y-
А формула C5) примет вид
N
Пусть N= 10. В качестве значений у используем тройки чисел
из таблицы А (умноженные на 0,001). Промежуточные результаты
сведены в табл. 1,
63

1 1
Y/
I)
singy
1
0,865
1,359
0,978
2
0,159
0,250
0,247
3
0,079
0,124
0,124
4
0,566
0,889
0,776
5
0,155
0,243
0,241
6
0,664
1,043
0,864
7
0,345
0,542
0,516
T
8
0,655
1,029
0,857
а б л i!
9
0,812
1,275
0,957
ца 1
10
0,332
0,521
0,498
Результат вычисления:
>	/ и 0,952.
' (б) Пусть теперь pg (х) =8 х/л2. Для разыгрывания | исполь-
используем уравнение B3)
г
Г (8лг/я2) dx — Y,
о
откуда после несложных вычислении получим
Формула C5) примет такой вид:
N
8N
Ij '
Пусть N=10. Числа у выберем те же, что в случае (а). Про-
Промежуточные результаты сведены в табл. 2.
Таблица 2
1
Y/
ll
sin Ij
ll
l
0,865
1,461
0,680
2
0,159
0,626
0,936
3
0,079
0,442
0,968
4
0,566
1,182
0,783
5
0,155
0,618
0,937
6
0,664
1,280
0,748
7
0,345
0,923
0,863
8
0,655
1,271
0,751
9
0,812
1,415
0,698
10
0,332
0,905
0,868
Результат вычисления:
/и 1,016.
Как мы и предполагали, второй способ вычисления дал более
точный результат,
54

8.4. Об оценке ошибки. В п. 8.2 уже отмеча-
отмечалось, что абсолютная ошибка при вычислении интегра-
интеграла / по формуле C5) практически не может превысить
величины 3]/DiiW. Однако в действительности ошибка,
как правило, оказывается заметно меньше этой величи-
величины. Поэтому для характеристики ошибки на практике
часто используют другую величину — так называемую
вероятную ошибку
бвер=0,675/бда
Фактическая абсолютная ошибка зависит от использо-
использованных в расчете случайных чисел и может оказаться
как в 2—3 раза больше, так	_.
и в несколько раз меньше,	ц
чем бвер- Так что бвер дает
нам не границу ошибки, а
порядок ее *).
Вернемся к примеру из п. 8.3.
По значениям, приведенным в
табл. 1 и 2, можно приближен-
приближенно сосчитать дисперсии Оц для
обоих методов расчета. Соответствующая расчетная формула при-
приведена в п. 6.1 **).
В табл. 3 для обоих методов расчета указаны приближенные
значения дисперсий Dt), сосчитанные по ним вероятные ошибки
бвер, а также фактические абсолютные ошибки, полученные при
расчете 6счет.
Мы видим, что 6Счет действительно имеют тот же порядок,
ЧТО бвер-
Метод
(а)
(б)
Dri
0,256
0,016
6вер
0,103
0,027
6счет
0,048
0,016
*) Определение вероятной ошибки см. в п. 9.7.
**) Для метода (а):
Ю	/ ю	\ 2-|
= |g- D,604 — 3,670) = 0,256.
Для метода (б):
я4
1 9 «Ь4
= =?g- F,875 — 6,777) = 0,016.
о/о

ДОБАВЛЕНИЕ
с
§ 9. Доказательства некоторых предложений
В этом параграфе приведены доказательства некото-
некоторых утверждений, высказанных в предыдущих парагра-
параграфах, которые кажутся нам
несколько громоздкими
для популярного изложе-
изложения или предполагают
больше знаний по теории
вероятностей.
9.1. Обоснование
метода Неймана
для разыгрывания
случайной величи-
величины (п. 4.3). Случайная
точка Г равномерно рас-
распределена в прямоуголь-
прямоугольнике abed (рис. 27), пло-.
щадь которого равна
М0(Ь—а)*). Вероятность
того, что точка Г окажет-
окажется под кривой у=р(х) и не будет отброшена, равна от-
отношению площадей
Р (х) dx
й)
*) Ср. п. 9.3.
56

А вероятность того, что точка Г окажется под кри-
кривой у=>р(х) на интервале а'<х<.Ь', также равна отно-
отношению площадей
у
J р (х) dx
Ма (Ь — а) Р
Следовательно, среди всех отобранных нами значе-
значений | доля значений, попавших в интервал (а1, Ь'), рав-
равна частному
Ь'
J p (x) dx	»•
Af, (b-a) : Мо (b - а) "" J P(x)dx>
a'
что и требовалось доказать.
9.2. Плотность распределения величины
t,' — a + at, (п. 4.4). Предполагается, что величина С нор-
нормальна с математическим ожиданием М? = 0 и диспер-
дисперсией D?=li так что ее плотность
рг (х) = -^=r e~J.
Чтобы вычислить плотность распределения величи*
ны Z,', выберем два произвольных числа Xi<.x2 и вычис-
вычислим вероятность
Р {*i < С < л-2} =Я {*, <
Следовательно,
(JTj-e)
В последнем интеграле сделаем замену переменной
х'=а + ах. Получим
2 (*'«)'
У 2л a J
57

откуда следует (ср. A4)) нормальность величины ?' с
параметрами М?' = а, D?' = or2.
9.3. Равномерное распределение точек
(y'> y") в квадрате (п. 4.5). Так как координаты точ-
точки (у', y") независимы, то плотность ее р(х,у) равна
произведению плотностей
р (х, у) = р ,{х) р „ (у).
Каждая из этих плотностей тождественно равна 1. Зна-
Значит, р($, у)^1 (при 0<х<1 и 0<у<1). А это и озна-
означает равномерность распределения точки (y', у") в еди-
единичном квадрате.
случайного направления
задавать направление единичным
вектором, выходящим из
начала координат. Концы
таких векторов располо-
расположены на поверхности еди-
единичной сферы. Слова «лю-
«любое направление одинако-
одинаково вероятно» означают,
что конец вектора пред-
представляет собой случай-
случайную точку Q, равномерно
распределенную на по-
поверхности этой сферы:
вероятность того, что п
окажется в любом эле-
эле9.4.
(п. 7,1).
др
Выбор
Условимся
Рис. 28.
dS,
менте поверхности
равна dS/4n.
Выберем на поверхности сферы сферические коорди-
координаты (ф, т|>) с полярной осью Ох (рис. 28). Тогда
dS = sin у d(fdi\	C6)
причем 0-<ф<я,
Так как координаты ф и -ф независимы, то плотность
точки (ф, 1|з) равна произведению /?(<р, гр) =рф (ф)/А1,-(^).
Из этого равенства, соотношения C6) и соотношения
/7(ф, -ф) dtp сГф
вытекает, что
58

Проинтегрируем это выражение по ф от 0 до 2л. Прини-
2я
мая во внимание условие нормировки | р^ (^) dty = 1,
о"
получим, что
C8)
Разделив C7) на C8) почленно, найдем, что
Очевидно, i[5 равномерно распределена в интервале
(О, 2л) и формула для разыгрывания i|) запишется так:
D0)
Формулу для разыгрывания ф найдем при помощи
уравнения B3):
ф
A/2) J
о
J
о
откуда
cos<p = l—2у.	D1)
Формулы D0) и D1) позволяют нам выбирать (ра-
(разыгрывать) случайное направление. Конечно, значения у
в этих формулах должны быть разными.
Формула D1) отличается от последней формулы
п. 7.1 только тем, что вместо у в ней стоит 1—у, а эти
величины имеют^одинаковое распределение.
9.5. Преимущество способа расчета с ве-
весами (п. 7.3). Введем случайные величины v и v', рав-
равные количеству (весу) прошедших сквозь пластинку ней-
нейтронов, полученному при расчете одной траектории ме-
методом п. 7.2 и методом п. 7.3.
По физическому смыслу задачи
Mv — Mv' = р+.
Так как v может принимать только два значения,
0 и 1, то распределение v задается таблицей .
1 0
59

Приняв во внимание, что v2 = v, нетрудно вычислить,
что Dv = p+— (p+f.
Легко видеть, что величина v' может принимать бес-
бесконечное число значений: wq=\, w\, w2, .. ., wh, ..., a
также значение 0. Поэтому ее распределение задается
таблицей
(w0 Wj w2 ... wh ... 0\
v =
Wo Я\ Яг ¦•¦ Як ¦•¦Я)
Значения q{ нас интересовать не будут, так как в любом
случае можно записать выражение для дисперсии
Заметив, что все Wh ^ 1 и что 2 1ЮкЯк== Mv' = p+,
получим неравенство Dv' < рь — (/>+)г = Dv.
Тот факт, что дисперсия v' всегда меньше, чем дис-
дисперсия v, показывает, что метод и. 7.3 для расчета р*
всегда лучше, чем метод п. 7.2.
Такой же вывод можно сделать о расчете /г, а при
не слишком большом поглощении — и о расчете р°.
9.6. Наилучший выбор! (п. 8.2). В п. 8.2 полу-
получено выражение для дисперсии Dt]. Чтобы найти мини-
минимум этого выражения при всевозможных выборах р\{х),
воспользуемся хорошо известным в анализе неравенством
Г»	-12	ft	b
\\\u(x)v{x)\dx\ < \a'1{x)dx- f v2(x)dx.
'	L a	-'a	a
Положим и = g{x)'YP> (x), а г) = У/?.(х); тогда из
этого неравенства получим
ь
Зна
g(x)]dx
чнт,
2 b
<f "'
а
(х
(х
- ь
j
- ft
b
a
\g(x)\dx
— P.	D2)
Остается доказать, что нижняя граница достигается
тогда, когда р\(х) пропорциональна |g(*)|,
60

Пусть
ь х) == „ lSK n		D3)
J | g (х) | dx
a
Нетрудно вычислить, что при такой плотности
ь
и дисперсия Dii действительно равна правой части D2).
Заметим, что выбрать «наилучшую» плотность D3)
для расчета практически невозможно: для этого нужно
знать значение интеграла j | g (x) | dx. А вычисление по-
а
следнего интеграла представляет собой задачу, равно-
равноценную стоящей перед нами задаче о вычислении инте-
ь
грала J g{x)dx. Поэтому мы ограничились рекоменда-
а
цией, указанной в п. 8.2.
9.7. Определение вероятной ошибки
Хп. 8.4). Пусть Z, — нормальная случайная величина, опре-
определенная в п. 2.3. Нетрудно вычислить, что, каковы бы
ни были а и а при г = 0,675 0,
а+г
f ft (jc) rf* = 0,5.
Отсюда вытекает, что
Я{|С — а\ <г}=Р{\1~а|>г}=-0.5,
то есть отклонение большее, чем г, и отклонение мень-
меньшее, чем г, одинаково вероятны. Поэтому величина г на-
называется вероятной ошибкой случайной величины ?.
В п. 8.1 мы вычисляем приближенно нормальную
величину р= A/N) (т]1 + т]2+ ... +г]Л-). Ее математиче-
математическое ожидание а=Мр = /, а дисперсия cr2=Dp= DnAV.
Поэтому вероятная ошибка величины р приближенно
равна 0,675 У[Ж

en
to
86515
69186
41686
86522
72 587
52 452
76 773
04 825
87113
84 754
0,2005
1,1609
0,5864
0,1425
0,9516
—0,5863
1,1572
—0,4428
—0,3924
0,8319
0,9780
90 795
03393
42 163
47 171
93 000
42 499
97 526
82134
84 778
57 616
1,1922
—0,6690
—0,9245
—0,2863
—1,7708
0,8574
0,9990
—0,5564
1,7981
0,4270
—0,7679
Таб
66155
42 502
85181
88 059
89 688
33 346
27 256
80 317
45 863
38132
Таб л
—0,0077
—1,5893
0,0904
1,2809
2,8854
—0,5557
—0,1032
—0,5098
0,6141
—0,8888
0,8960
ТА
БЛИЦЫ
лица А. 400 случайных
66 434
99 224
38 967
89 342
78 416
83 935
66 447
75 120
24520
64 294
ица Б. 88
0,0348
0,5816
1,5068
0,4043
0,4686
0,8115
0,5405
—1,1929
—1,3596
0,4167
0,5154
*) Случайные цифры имитируют значения случайной !
••) Нормальные величинь
а — О, а = 1.
имитируют
56 558
88 955
33181
67 248
27 589
79 130
25 731
45 904
19 976
15 218
Цифр *)
12 332
53 758
72 664
09 082
99 528
90 410
37 525
75 601
04 925
49 286
нормальных величин **)
1,0423
1,8818
—1,1147
0,6379
1,4664
—0,2676
—0,6022
—0,0572
1,4943
—0,8513
—0,7165
—1,8149
0,7390
0,2776
—0,4428
1,6852
—1,2496
0,0093
—0,5061
—0,4406
1,1054
0,8563
94 377
91641
53807
12 311
14 480
45 420
16 287
70 492
07 824
89571
1,1803
—0,2736
0,1012
—2,3006
—0,9690
—1,2125
0,2119
—0,1557
—0,2033
1,2237
—1,1630
!еличины с распределением B2) (см. п. 3.1).
значения нормальной (гауссо
вской) случайной величины
57 802
18 867
00 607
90316
50 961
77 757
66181
10 274
76 044
42 903
0,0033
1,0828
—1,3566
—0,6446
—0,0831
1,3846
—1,4647
—1,2384
—0,1316
—0,7003
1,8800
? с параметрами

ЛИТЕРАТУРНЫЕ УКАЗАНИЯ
Наиболее полное изложение вопросов, связанных с методом
Монте-Карло, имеется в книгах:
[1] Бусленко Н. П., Ш рей дер Ю. А., Метод статистических
испытаний (Монте-Карло) и его реализация в цифровых маши-
машинах, Физматгиз, 1961.
[2] Бусленко Н. П., Голенко Д. И., Соболь И. М., С р а-
г о в и ч В. Г., Ш р е й д е р Ю. А., Метод статистических испыта-
испытаний (метод Монте-Карло), Физматгиз, СМБ, 1962.
Статьи по многим вопросам, из числа затронутых в настоящей
книжке, имеются в энциклопедии:
[3] Автоматизация производства и промышленная электроника,
тт. 1—4, БСЭ, 1962—1965.
Можно указать также некоторые книги, посвященные отдель-
отдельным вопросам:
К § 4
[4] Г о л е н к о Д. И., Моделирование и статистический анализ псев-
псевдослучайных чисел на ЭВМ, «Наука», 1965.
К § 5
[5] Б у с л е н к о Н. П., Математическое моделирование производ-
производственных процессов на цифровых вычислительных машинах,
«Наука», 1964.
К § 6
1етоды
ргия», 1
По вопросам, относящимся к §§ 7, 8, см. книгу [2].
[6] Т у р ке л ьт а у б Р. М., Методы исследования точности и на-
надежности аппаратуры, «Энергия», 1966.

ОГЛАВЛЕНИЕ
Предисловие		3
Введение	 . • •	5
§ 1. Общее представление о методе		5
Глава 1. Моделирование случайных величин ........	9
§ 2. Случайные величины		9
§ 3. Получение случайных величин на ЭВМ		20
§ 4. Преобразования случайных величии		25
Глава 2. Примеры применения метода Моите-Карло ....	33
§ 5. Расчет системы массового обслуживания		33
§ 6. Расчет качества и надежности изделий		38
§ 7. Расчет прохождения нейтронов сквозь пластинку . .	44
§ 8. Вычисление определенного интеграла		51
Добавление		56
§ 9. Доказательства некоторых предложений		56
Таблицы		62
Литературные указания		63

Цена 10 коп.
Вьш
Вып.
Вьш
Вып.
Вьш
Вып.
Вьш.
Вып.
Вып.
Вып.
Вып.
Вьш.
Вып.
Вып.
Вьш.
Вьш.
Вып.
Вьш.
Вып.
Вып.
Вып.
Вып.
Вьш.
Вьш.
Вьш.
Вып.
Вьш.
Вып.
Вьш.
Вып.
Вьш.
Вып.
Вып.
Вьш.
Вып.
Вып.
Вып.
Вьш.
Вып.
Вып.
Вып.
Вып.
Вып.
Вьш.
Вып.
1
2.
3
4.
5
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
ПОПУЛЯРНЫЕ ЛЕКЦИИ ПО МАТЕМАТИКЕ
А.
И
И
А.
П
Н
А.
и
п.
с
и.
п
и
г.
пеней
А.
А.
А.
Я.
И.
А.
о.
и.
с.
с.
п
и
Маркушевич. Возвратные последовательности.
Натансон. Простейшие задачи на максимум и минимум.
Соминский. Метод математической индукции.
Маркушевнч. Замечательные кривые.
Коровник. Неравенства.
Воробьев. Числа Фибоначчи.
Курош. Алгебраические уравнения произвольных сте-
Гельфонд. Решение уравнений в целых числах.
Маркушевич. Площади н логарифмы.
Смогоржевский. Метод координат.
Дубнов. Ошибки в геометрических доказательствах.
Натансон. Суммирование бесконечно малых величин.
Маркушевич. Комплексные числа л_-КОН_фовмные_ото-
браження. *~
А. И. Фетисов. О доказательствах в геометрии.
И.
В.
В.
г.
л.
А.
Л.
В.
А.
Б.
А.
Б.
В.
р.
г.
г.
м.
А.
м.
и.
г.
с.
и.
с.
А.
А.
Шафаревич. О решении уравнений высших степеней.
Шерватов. Гиперболические функции.
Болтянский. Что такое дифференцирование?
Мнракьян. Прямой круговой цилиндр.
Люстерник. Кратчайшие линии.
Лопшиц. Вычисление площадей ориентированных фигур.
Головина и И. М. Яглом. Индукция в геометрии.
Болтянский. Равновеликие н равносоставленные фигуры.
Смогоржевский. О геометрии Лобачевского.
Аргунов н Л. А. Скорняков. Конфигурационные теоремы.
Смогоржевский. Линейка в геометрических построениях.
Трахтенброт. Алгоритмы н машинное решение зад^ч.
Успенский. Некоторый прГО1бж"ёТШ[—мтагннящ-ТГмате-
матнке
Н.
А.
Архангельский и Б. И. Зайцев. Автоматические цнфро-
А.
н.
КОстовский. Геометрические построения одним цнр-
кулем.
Г.
А.
Е.
А.
Б.
Н.
В.
Е.
Г.
С.
С.
Е.
Я.
Г.
Шилов. Как строить графики.
Дорфман. Оптика конических сечений.
Вентцель. Элементы теорнн игр.
КЯПГОВ. ЧтП ТЯКПР J|H4fUHnfi' ПрпгрячМ1фП1>яиио
Маргулис. Системы линейных уравнений.
Виленкин. Метод последовательных_щ}нблнженнй.
Болтянский. "Огибающая.
Г. Е. Шилов. Простая гамма (устройство музыкальной
Ю.
Н.
С.
Б.
Ю.
А.
Н.
В.
Ю.
И
Шрейдер. Что такое расстояние?
Воробьев. Признаки делимости.
Фомин. Системы счисления.
Коган. Приложение механики к геометрии
Любич н Л. А. Шор. Кинематический метод в геомет-
рнческнх задачах.
В.
И.
И.
А.
Я.
М.
Успенский. Треугольник Паскаля.
Бакельман. Инверсия.
Яглом. Необыкновенная алгебра.