Материал: Имитационное и статистическое моделирование — копия
3.2. Моделирование непрерывных случайных величин
Для построения датчиков непрерывных с. в., имеющих заданное распределение вероятностей, разработано большое число алгоритмов [2, 6–15]. В современной учебной литературе по моделированию методы построения датчиков непрерывных с. в. наиболее полно разбираются в книге [8]. Основные методы реализации непрерывных с. в. разъясняются также в [16]. Вместе с тем, во многих языках и пакетах моделирования имеется большое число удобных для использования встроенных датчиков непрерывных с. в.
С учётом сказанного ниже приводятся лишь алгоритмы и формулы для реализации наиболее часто встречающихся непрерывных с. в. и некоторые примеры их реализации на Excel.
3.2.1. Моделирование равномерной случайной величины
С. в. x ~ R[A, B], равномерно распределённую на интервале A ≤ t ≤ B, можно реализовать путём линейного преобразования БСВ z:
x = (B – A)·z + A. |
(1.12) |
При этом из (1.12) получаем |
|
M(x) = A + (B – A)/2, |
(1.13) |
D(x) = (B – A)2/12. |
(1.14) |
3.2.2. Моделирование экспоненциальной случайной величины |
|
Экспоненциальная с. в. x имеет ф. р. в. |
|
F(t) = 1 – e−λt , |
(1.15) |
где t ≥ 0, параметр λ > 0 . Её м. о. и дисперсия |
|
M(x) = 1/λ , D(x) = 1/λ2. |
(1.16) |
Экспоненциальную с. в. можно реализовать с помощью следующего пре- |
|
образования БСВ z: |
|
x = – (1/λ)·ln(z). |
(1.17) |
Пример 7 . На чистом листе Ms Excel в диапазоне A1:A1000 с помощью функции =СЛЧИС() создайте выборку из 1000 значений БСВ. Полагая λ = 1, в соответствии с (1.17) запишите в ячейке B1 формулу =–Ln(A1). Двойным щелчком по маркеру заполнения скопируйте её в диапазон B1:B1000. В этом диапазоне получена выборка с. в. x с экспоненциальным распределением (1.15) при λ = 1.
Проверьте вид распределения полученной с. в. путём построения её эмпирической ф. р. в. Для этого скопируйте диапазон B1:B1000 в буфер и вставьте в C1 посредством меню ...Специальная вставка/значения, после чего отсортируйте диапазон C1:C1000 по возрастанию. В ячейки диапазона D1:D1000 запишите арифметическую прогрессию 0,001, 0,002, ... , 1,000 и постройте диаграмму Точечная по диапазону C1:D1000, как в п. 2.4. Диаграмма будет выглядеть примерно как на рис. 1.3 («дрожащая» кривая линия, фрагмент которой справа увеличен).
16
Рис. 1.3. Сравнение эмпирической и теоретической ф. р. в. экспоненциальной с. в.
Гладкую теоретическую кривую, приведенную на рис. 1.3, можно добавить к вашей диаграмме, если в ячейку E1 ввести формулу (1.15) для вычисления теоретического значения ф. р. в. в виде =1-exp(-C1); затем её надо скопировать на весь диапазон E1:E1000. Добавление этого диапазона данных к диаграмме производится через контекстное меню диаграммы Исходные данные.../Ряд/
Добавить.
3.2.3. Моделирование эрланговской случайной величины
Эрланговская с.в. x порядка k ≥ 1 имеет ф. р. в.
|
|
k −1 |
|
i |
|
|
|
||
F(t) |
= 1−∑ |
(λt) |
|
e−λt |
, |
t ≥ 0, |
(1.18) |
||
|
|
||||||||
|
|
i=0 |
i! |
|
|
|
|
|
|
и п. р. в. |
|
λ(λt)k −1 |
|
|
|
|
|
|
|
f (t) = |
e−λt , |
|
t ≥ 0, |
(1.19) |
|||||
(k −1)! |
|
||||||||
|
|
|
|
|
|
|
|
||
где параметр λ > 0 . Её м. о., дисперсия и начальные моменты r-го порядка таковы:
M(x) = k/λ , |
D(x) = k/λ2, |
M(xr) = k(k+1) ··· (k+r–1)/λr. |
(1.20) |
Поскольку распределением Эрланга обладает сумма k независимых экспоненциальных с. в., имеющих одно и то же значение параметра λ, то, с учётом (1.17), сгенерировать эрланговскую с. в. x можно просто как сумму:
k |
|
1 |
|
|
1 |
|
|
|
|
x = ∑i=1 |
− |
ln zi |
= − |
ln(z1 |
... zk ) , |
(1.21) |
|||
λ |
λ |
где zi (i = 1, ..., k) – независимые реализации БСВ.
3.2.4. Моделирование нормальной случайной величины
Утверждение, что некоторая с. в. x имеет нормальное распределение с м. о. M(x) = μ и дисперсией D(x) = σ2, записывают в виде x ~ N(μ, σ). П. р. в.
этой с. в.
17
f (t )= |
1 |
e− |
(t −μ )2 |
|
|
2σ 2 . |
(1.22) |
||||
2πσ |
|||||
|
|
|
|
Её ф. р. в. F(t) не выражается в элементарных функциях в виде конечной формулы. Для реализации любой нормальной с. в. достаточно иметь датчик стандартной (т. е. нормированной и центрированной) нормальной с. в. x ~ N(0,1) . Чтобы реализовать с. в. x с распределением (1.22) используют следующее линейное преобразование стандартной нормальной с. в. x :
x = σ x + μ. |
(1.23) |
При этом стандартную нормальную с. в. часто реализуют приближённо, как сумму других с. в., основываясь на центральной предельной теореме теории вероятностей. Например, её можно реализовать в виде такой суммы двенадцати независимых значений БСВ:
12
x = ∑zi − 6. (1.24)
i=1
Однако такой подход даёт плохое приближение для больших уклонений от среднего, превышающих 2σ.
Метод Бокса и Мюллера позволяет получить два независимых значения x1 и x2 стандартной нормальной с. в. из двух независимых значений z1 и z2 БСВ по формулам:
x |
|
= |
− 2 ln z |
1 |
sin(2π z |
2 |
), |
|
|||
|
1 |
|
|
|
|
|
(1.25) |
||||
|
|
|
= |
− 2 ln z |
|
cos(2π z |
|
). |
|||
x |
2 |
1 |
2 |
|
|||||||
|
|
|
|
|
|
|
|
|
|||
Этот метод точный, но считается трудоёмким [8]. Однако, как показывают эксперименты, это мнение устарело ввиду того, что современные персональные компьютеры оснащены арифметическими сопроцессорами. Точный метод (1.25) в действительности оказывается также и более быстрым, чем приближённый метод (1.24).
Приведём ещё формулу для чётных начальных моментов r-го порядка с. в. x :
M (x 2k ) =(2k–1)!! = (2k–1)·(2k–3)· … ·1. |
(1.26) |
Все её моменты нечётного порядка равны нулю. 3.2.5. Применение гистограмм и полигонов
В [22] приводятся алгоритмы реализации тех с. в., у которых п. р. в. имеют форму гистограммы (рис. 1.4, а) или полигона (рис. 1.4, б). Такие п. р. в. появляются обычно либо как результат обработки эмпирических данных, либо как результат аппроксимации теоретических распределений. Общая площадь, заключённая внутри гистограммы или полигона, равна единице (свойство п. р. в.).
Реализацию с. в. по заданной гистограмме можно осуществить путём случайного выбора номера интервала постоянства п. р. в. (вероятность каждого такого интервала гистограммы равна площади построенного на нём прямоугольника) и последующей реализации на этом интервале равномерной с. в.
18
f(t) |
f(t) |
а) |
б) |
t |
t |
Рис. 1.4. Гистограмма (а) и полигон (б) |
|
Если п. р. в. f(t) задана в виде полигона, то можно перейти от неё к ф. р. в. F(t) и разыграть с. в. x методом обращения [16]. Поскольку F(t) является интегралом от п. р. в.:
F (t) = ∫t |
f (s)ds , |
(1.27) |
−∞ |
|
|
то здесь F(t) представляет собой кусочно-квадратичную зависимость, и её график агрегируется из отрезков разных параболических дуг. Тем не менее из непрерывности f(t) следует, что график функции F(t) – гладкий, без изломов. Построенный методом обращения этой ф. р. в. генератор с. в. x должен содержать процедуру, которая по разыгранному значению БСВ z определяет на оси ординат графика F(t) тот интервал значений, на который проецируется одна из агрегированных дуг и в который попало значение z. После этого x можно вычислить как значение обратной от z функции
x = F −1 (z) |
(1.28) |
по квадратному уравнению этого отрезка дуги.
Пример 8 . Предположим, график f(t) представляет собой две боковые стороны равнобедренного треугольника с основанием на отрезке (0, 2). Из равенства площади треугольника единице следует, что его высота равна 1. Из (1.27) по этой f(t) определяем:
0,
t 2 F (t) = 2t1,
|
t < 0, |
|
/ 2, |
0 ≤ t <1, |
(1.29) |
−t 2 / 2 −1, |
1 ≤ t < 2, |
2 ≤ t.
Применение к этой ф. р. в. обращения (1.28) приводит к следующей комбинированной формуле для генерации x:
|
2z , |
z <1/ 2, |
|
|
(1.30) |
||
x = |
|
|
|
|
− 2 − 2z , |
1/ 2 ≤ z. |
|
2 |
|
В связи с рассмотренным примером заметим, что описанную в нём «треугольную» п. р. в. имеет также сумма двух независимых БСВ. Следовательно, вместо (1.30) для генерации x можно использовать более простую формулу x = z1 + z2.
19
Это наблюдение позволяет предположить, что для генерации с. в., которые задаются кусочно-линейными п. р. в. – полигонами, можно найти метод, более эффективный, чем метод обращения, если использовать суммирование равномерных или кусочно-равномерных случайных величин. По всей видимости, наряду с операцией суммирования в этом методе целесообразно было бы использовать операцию случайного выбора (т. н. смешивание) случайных величин.
Проверку этого предположения оставляем читателю.
3.3. Моделирование случайных процессов
Всякий случайный процесс – дискретный или непрерывный – моделируют последовательностью с. в., которые генерируются либо как статистически независимые, либо как взаимозависимые величины. При этом применяют разные способы реализации зависимости между с. в. Чаще всего, поскольку модель всегда имеет какой-либо «физический смысл», способ реализации зависимости между с. в. определяется причинно-следственными связями, действующими в моделируемой системе. Во многих случаях идеализация и упрощение моделируемых причинно-следственных отношений приводит к типичным математическим схемам описания случайных процессов. К числу таких типичных математических схем относятся Марковские процессы, Гауссовы процессы, Броуновское движение, случайное блуждание на прямой и др. [23, 24].
Свойства процессов, описываемых подобными схемами, хорошо изучены методами теории вероятностей. Но в тех случаях, когда в отношении этих процессов решаются задачи, отличающиеся от классических задач, или когда сами исследуемые процессы не вписываются в изученные математические схемы, используют статистическое моделирование. В качестве примера такого статистического моделирования рассмотрим моделирование Броуновского движения.
3.3.1. Броуновское движение на неограниченной плоскости
Процесс Броуновского движения на плоскости можно приближённо представить в виде последовательности случайных перемещений точечной частицы, выполняемых в дискретные моменты времени t = 0, 1, 2, … .
Будем полагать, что в момент t = 0 частица находится в начале системы координат. Перемещение частицы в любой из следующих моментов времени реализуется выбором случайного направления в пределах от 0 до 2π радиан и смещением частицы в этом направлении на шаг единичной длины.
Модель такого Броуновского движения легко можно реализовать средст-
вами Ms Excel.
Пример 9 . Откройте чистый лист Excel, введите в ячейки A1, B1 и C1 три заголовка Угол, X, Y и под ними (соответственно в ячейках A2:C3) введите формулы:
Угол |
X |
Y |
=2*ПИ()*СЛЧИС() |
=COS(A2) |
=SIN(A2) |
=2*ПИ()*СЛЧИС() |
=B2+COS(A3) |
=C2+SIN(A3) |
Нетрудно видеть, что непосредственно под заголовками полученной таблички вычисляется положение частицы в момент времени t = 0. В первой ко-
20