CC BY
85
№ 5 (35) 2011
А. А. Емельянов, докт. экон. наук, профессор МФПУ «Синергия»
Лаг-генераторы для моделирования рисковых ситуаций в системе Actor Pilgrim
Создание программ имитационного моделирования ведется, как правило, с помощью языков моделирования и моделирующих комплексов. Важным компонентом последних являются генераторы случайных чисел и последовательностей, требования к которым делают решение задачи их создания далеко не тривиальным.
Введение
Имитационное моделирование довольно успешно используется при исследовании рисковых ситуаций [1, 3]. Однако компьютерные имитационные модели могут привести к неточному результату, если исследователи на этапе их создания не уделят должного внимания двум важным вопросам обеспечения точности моделирования:
1) как спланировать времяемкий компьютерный эксперимент:
• чтобы с высокой точностью получать значения и координаты экстремальных точек при исследовании рискового процесса;
• избежать клинча таймера имитационной модели, когда элементарные отрезки времени в процессе выполнения модели становятся пренебрежительно малыми по сравнению с накопленным значением таймера и соизмеримыми с погрешностями численных вычислений и поэтому перестают «правильно двигать» модельное время при регистрации всех моделируемых событий, включая как редкие рисковые события, так и регулярные, часто происходящие;
2) какими программными средствами обеспечить получение корректных последовательностей псевдослучайных величин.
На первый вопрос можно найти подходящие ответы или способы самостоятельного получения решений в научной и учебно-методической литературе [1, 3, 5]. Но на вто-
рой вопрос практически ни в одном учебнике готовых ответов нет. Дело в том, что программные генераторы получения последовательностей псевдослучайных чисел давно используются на практике, в последнее время — в имитационном моделировании рисковых процессов [3]. Причем почти все используемые стандартные библиотечные программы-генераторы основаны на конгруэнтных методах, разработанных Д. Лемером1 [8]. Эти программы характеризуются:
• значительным быстродействием и компактностью в смысле требований к памяти компьютера;
• величиной периода получаемой последовательности в пределах от 1/2 до 1 миллиарда чисел.
Как правило, конгруэнтные программы-генераторы оперируют с 32-разрядными машинными словами. Свойства этих генераторов подходят для большинства научных и учебных применений. Однако особенности распределения случайных чисел, генерируемых линейным конгруэнтным алгоритмом (главным образом — период псевдослучайной последовательности и неравнозначность битов получаемого случайного числа по статистическим свойствам), в отдельных случаях затрудняют и даже делают невозможным использование таких генераторов
1 Lehmer D. H. Mathematical Methods in Large-Scale Computing Units. Annals Computer Laboratory, Harvard University. V. 26. 1951.
№ 5 (35) 2011
в имитационных моделирующих программах, требующих высокого разрешения.
Типичными примерами задач, где требуется повышенная точность генераторов случайных числовых последовательностей, являются следующие:
• исследование редких рисковых событий на длинной оси времени при одновременном моделировании хозяйственно-экономической, финансовой, военной или иной деятельности, для которой эти события представляют угрозу;
• изучение длинных «хвостов» распределений вероятностей рисковых событий для получения возможностей составления прогнозов.
В таких задачах наряду с отлаженными в течение многих лет модификациями конгруэнтных алгоритмов хорошие перспективы в настоящее время имеются для семейства фибоначчиевых алгоритмов [10]. Их практическое применение стало возможным после появления 64-разрядных процессоров, обладающих высокой скоростью выполнения операций вещественной арифметики. Несмотря на сложность их программной реализации, такие алгоритмы можно рекомендовать в качестве инструментария проведения исследований,критичных к качеству случайных чисел.
Мультипликативный конгруэнтный метод
Конгруэнтный метод представляет собой арифметическую процедуру для генерирования конечной последовательности равномерно распределенных чисел. С использованием нескольких рекуррентных формул было построено множество конгруэнтных алгоритмов. В основе каждого из них лежит фундаментальное понятие конгруэнтности [8, 9, 15].
Определение. Два натуральных числа А и В конгруэнтны (сравнимы) по модулю т, где т — целое число, в случае, если существует такое целое число к, что А - В = к • т, то есть если разность А - В делится на т,
и числа А и В дают одинаковые остатки от деления на m. Это определение записывается как
A = B mod m
и читается «А конгруэнтно B по модулю m».
Наиболее известны три разновидности конгруэнтных алгоритмов: мультипликативный, смешанный и аддитивный. В большинстве стандартных программ используется мультипликативный алгоритм. Со смешанным и аддитивным методами можно подробно ознакомиться в [8].
Итеративная формула мультипликативного конгруэнтного метода имеет вид:
xi = (cxi_1 )mod m ,
(1)
где c — положительное целое число.
Известны и другие разновидности формулы (1) для реализации конгруэнтного метода [8, 15].
Сущность мультипликативного метода заключается в многократном выполнении операции умножения на специально подобранную константу c (в связи с чем он и получил свое название). Для получения нового псевдослучайного числа через последнее значение xj мы должны сделать изменение типа j = j + 1. После этого существующее число обозначается как х _ 1. Далее нужно взять последнее случайное число х1 - 1, умножить его на постоянный коэффициент с и взять модуль полученного числа по т. Это означает, что произведение с • xj- 1 делится на т, и остаток считается новым псевдослучайным числом х. Поэтому для генерирования последовательности чисел х1 необходимо знать множитель с и модуль т, а также начальное значение последовательности х0 которое на первом шаге работы алгоритма будет использоваться в качестве последнего значения х.
Любой генератор псевдослучайных чисел может дать лишь конечное множество целых случайных чисел, после которого последовательность начинает повторяться.
I
Uj
99
-ч ПРИКЛАДНАЯ ИНФОРМАТИКА
№ 5 (35) 2011 ' -
.«а <Е
S
и ^
i
is
IS
an 1
со
Si
и
I
S 1
i
i
is §
Период или длина последовательности T зависит от разрядной сетки вычислительной машины и от величины выбранного модуля m, а статистические свойства полученной последовательности зависят от выбора начального значения x0 и множителя с. Таким образом, выбирать с, x0 и m следует так, чтобы обеспечить максимальный период и минимальную корреляцию между генерируемыми числами.
Правильный выбор модуля не зависит от системы счисления, используемой в данном компьютере. Самым естественным было бы взять m равным длине машинного слова и работать с целыми числами. При использовании компьютера, где применяется двоичная система счисления, m будет равно 2', где l — число разрядов (бит) в машинном слове. Максимальный период при правильном (удачном) выборе с и x0 определяется соотношением T = 2'- 2 или T = m/4. Период последовательности для 32-разрядного целого числа равен T = 1073741824, т. е. порядка одного миллиарда. Но если один разряд из 32 расходуется на знак числа, т. е. не используется в процедуре генерации, то максимальное значение периода равно T = 536870912, т. е. примерно полмиллиарда чисел. Такая длина может быть вполне удовлетворительной для решения многих учебных и исследовательских задач. Однако для ряда специальных приложений, связанных с численными исследованиями в области теории рисков, этого недостаточно.
Количество модификаций конгруэнтных генераторов сейчас довольно велико [8]. Причины модификаций — борьба либо за увеличение периода, либо за улучшение статистических характеристик генератора. Наиболее удачным с точки зрения статистических характеристик подобия равномерному распределению можно считать 32-разрядный генератор Р. Шеннона [15]. Тесты показывают, что он, как по своим статистическим характеристикам, так и по длине периода существенно превосходит программный генератор, используемый в табличном процессоре Excel.
Пример программной реализации генератора приведен ниже.
long x; ке,
// Число на 3 2-разрядной сет-
double shannon(void) {
long y; // Рабочая переменная
double randum; // Число на интервале [0,1]
y = x * 1220703125; if (y < 0)
y = y + 2147483647 + 1; randum = y * (double)0.46566128730773 92578125e-9; x = y;
return randum; }
Очередное псевдослучайное число получается функцией shannon в переменной long x после каждого нового обращения к ней.
Данная функция нуждается в следующих пояснениях:
1) глобальная 32-разрядная переменная x типа long используется для получения псевдослучайных чисел x,, равномерно распределенных по всему диапазону чисел, определяемой разрядной сеткой, после каждого обращения к функции shannon;
2) величина множителя с = 1220703125 = = 513 была подобрана Р. Шенноном;
3) величина рабочей переменной увеличивается на 2147483647 + 1 в том случае, если возникает переполнение разрядной сетки после очередного умножения на с (десятичное число 2147483647 обеспечивает заполнение единицами всех разрядов целого двоичного числа);
4) в переменной randum с плавающей точкой типа double получается нормированное псевдослучайное число x, распределенное на интервале [0,1].
Переход от целого числа, распределенного на 32-разрядной сетке, к числу с плавающей точкой, распределенному на интервале [0,1], с учетом знакового разряда осуществляется умножением на константу 2- 31:
(double)0.4656612873077392578125e-9.
№ 5 (35) 2011
В качестве начального значения x0 для этого датчика рекомендуется выбрать число, которое не превосходит 109 и не делится на 2 или 5 (рекомендовано в [15]). При отладке моделирующих программ, использующих датчик, для получения конкретной цифровой последовательности x0 можно использовать таблицы случайных величин или, например, цифровую последовательность 123456789, которая удовлетворяет приведенной выше рекомендации.
При работе с отлаженной программой можно «усилить» иллюзию случайности последовательности, используя в качестве x0 кодовую последовательность таймера компьютера, который всегда имеет разные значения:
x0 = (long)time(NULL);
while ( (x0 % 2 == 0) || (x0 % 5 == 0) )
x0++; .
Следует отметить, что функция time в современных компьютерах возвращает целое число, которое представляет собой величину интервала времени (в секундах), начавшегося в 00 часов 01 января 1970 года и закончившегося моментом выполнения функции time в компьютере. Этот интервал представляет собой 64-разрядное число (_int64). Чтобы 32-разрядная переменная x0 включала все разряды таймера, перед time записывается определитель типа long. Понятно, что практически после каждого нового запуска моделирующей программы величина x0 будет получать новое, ранее не использованное значение. В результате и последовательности чисел, получаемых при разных запусках, оказываются различными.
Описанный выше конгруэнтный датчик обладает высоким быстродействием и лучшими статистическими характеристиками в семействе 32-разрядных конгруэнтных датчиков. Значение множителя c = 1220703125 и правила выбора начального значения x0 гарантируют то, что длина периода T будет не меньше, чем полмиллиар-
да. Однако и этой величины периода в не- | которых случаях недостаточно. г§
I
Методы увеличения периода ^
псевдослучайной последовательности
В настоящее время есть два основных направления увеличения периода:
1) увеличение модуля m конгруэнтного датчика;
2) установление зависимости очередного случайного числа x, сразу от нескольких параметров, а не только от его предшествующего значения как в конгруэнтных датчиках.
Увеличение модуля конгруэнтного датчика. Этот прием можно выполнить программным способом, искусственно получая длинное слово из нескольких машинных слов. Другого способа, особенно для 32-х разрядных компьютеров, нет. Однако этот прием приводит к резкому ухудшению быстродействия программного датчика, что может отрицательным образом сказаться в случаях больших (порядка многих миллиардов) значениях периода T. Поскольку современные процессоры в компьютерах работают с 64-х разрядными, а в суперЭВМ со 128-разрядными словами, на первый взгляд, напрашивается простое решение: увеличить период T за счет формального перевода программного датчика shannon на работу, например с 64-разрядными целыми числами (_int64). Для этого нужно только переопределить переменные x, y и преобразовать константу-множитель в константу типа double:
_int6 4 x;
_int6 4 y;
(double) 1.08420217248550443400745280086 9 9e-19; ,
где константа-множитель представляет собой величину 2 - 63.
Но, к сожалению, подобрать величины множителя c и начального значения x0, обеспечивающие такие же хорошие статистиче-
л 101
№ 5 (35) 2011
ские характеристики 64-разрядного генератора, какими обладает 32-разрядная версия программного генератора shannon, не удается: они намного хуже. Существует семейство тестов на проверку случайности в получаемой последовательности чисел, с которыми можно ознакомиться у Д. Кнута [9] и Р. Шеннона [15]. В соответствии с некоторыми из них в 64-разрядном датчике, который мы только что построили, имеет место существенная неравнозначность битов получаемого случайного числа по статистическим свойствам, т. е. увеличение периода T привело к тому, что получаемая последовательность чисел стала далека от равномерного распределения. Поэтому именно описанный выше 32-разрядный генератор Р. Шеннона пользуется популярностью в большинстве типовых программных приложений.
Определение очередного числа через два предшествующих значения. Простейшей последовательностью, в которой x, зависит более чем от одного из предыдущих значений, является последовательность чи-
== сел Фибоначчи о.
■§ X: = (x,_1 + x, 2) modm . (2)
! В отличие от конгруэнтного датчика, ис-<3 пользующего формулу типа (1), здесь ис-^ пользуется комбинация из двух предшест-Ц вующих чисел псевдослучайной последо-Ц вательности. Данный генератор рассматри-5 вался в начале 50-х гг., и обычно он дает | длину периода, большую, чем т. Тесты, од-| нако, показывают, что числа, получаемые ¡1 с помощью рекуррентного соотношения Фи-| боначчи недостаточно случайны. Таким об-| разом, выражение (2) можно рассматривать ss как элегантный, но неудачный пример источника случайных чисел. § Намного лучший аддитивный генератор был изобретен Дж. Ж. Митчелом и Д. Ф. Му-|| ром [8]. Они предложили несколько не-| обычную последовательность, определении ную так:
¡в x, = (x,_24 + x,_ss) modm , (3)
где i > 55, т — четное число, а x0, ..., x54 — произвольные целые не все четные числа. Числа-смещения 24 и 55 в данном определении были получены авторами для обеспечения периода последовательности xt не хуже, чем 255 - 1. В [9] доказывается, что при работе с целыми числами и /-разрядными машинными словами, когда модуль m = 2/, длина периода последовательности, определенной в выражении (3), точно равна
T = (255 _ 1)- 2/_1,
т. е. намного больше, чем у конгруэнтного генератора.
Применение алгоритма с лагами запаздывания
В отличие от рассмотренной выше последовательности чисел Фибоначчи алгоритм Фибоначчи с лагами запаздывания дает возможность получать очень длинные псевдослучайные последовательности, обладающие высокими статистическим характеристиками. Фибоначчиевы генераторы на базе этого алгоритма реализуются в вещественной арифметике. Перспектива использования таких генераторов появилась сравнительно недавно [10] в связи с совершенствованием материальной базы информатики (hardware): в процессорах современных компьютеров скорость выполнения арифметических операций с вещественными числами сравнялась со скоростью целочисленной арифметики. Это относится как к числам float (с плавающей точкой и одинарной точностью, длиной 32 бита), так и к числам double (с плавающей точкой и двойной точностью, длиной 64 бита). Далее будем иметь дело с числами типа double. Однако при реализации такого генератора ничто не мешает вместо вещественных чисел из диапазона [0,1] брать 32-разрядное целое число.
Число double с двойной точностью. Числа двойной точности с плавающей точкой занимают в памяти компьютера 64 бита. Как правило, в процессорах современ-
ПРИКЛАДНАЯ ИНФОРМАТИКА /-
' № 5 (35) 2011
ных компьютеров формат таких чисел соответствует стандарту IEEE 754 - 2008 [13] (рис. 1).
Экспонента (11 бит)
Мантисса (52 бита)
63 56i55
!i47 40i39 3231
йгРЧ!
115 8i7
Рис. 1. Формат числа с плавающей точкой double
Такие числа имеют точность около 16 десятичных значащих цифр и порядки в диапазоне примерно от 10 -308 до 10308. В компьютерах, которые работают с 64-разрядными числами с плавающей точкой, большинство численных вычислений осуществляется с двойной точностью, поскольку использование 32-разрядных чисел одинарной точностью (типа float) практически не повышает производительность компьютера, но ухудшает точность расчетов.
Вычислительная рекуррентная формула. Один из фибоначчиевых генераторов основан на следующей формуле [7]:
если xt a > x, b
x =i - ' " i-a >i-b , (4)
IXi-a-Xi-b + 1 еСЛИ Xi-a > Xi-b
где Xj — вещественные числа из диапазона [0,1];
a, b — целые положительные числа, называемые лагами запаздывания; i = max {a, b} + 1, max {a, b} + 2, max {a, b} + 3______
Для практической реализации фибонач-чиевого датчика на основе целых 32-разрядных чисел достаточно формулы
но при этом нужно учитывать эффект арифметических переполнений.
Наличие двух лагов запаздывания в вычислительной формуле стало причиной, по которой генераторы этого класса получили название лаг-генераторов (англ. laggenerators или lagged Fibonacci generators).
В англоязычной литературе фибоначчиевы генераторы иногда называют SWBG (sub-tract-with-borrow generators [10]).
При реализации фибоначчиевого генератора через целые 32-разрядные числа достаточно формулы
Xi = Xi-a ~Xi-b,
но при этом нужно учесть, что будут происходить арифметические переполнения, с которыми нужно что-то делать.
Получаемые таким генератором случайные числа обладают хорошими статистическими свойствами, причем все биты случайного числа равнозначны в статистическом смысле.
Значения лагов. Лаги a и b — «магические», и их не следует выбирать произвольно [7, 9]. Рекомендуются следующие значения пары лагов (a, b):
(17, 5), (20, 5), (55, 24) и (97, 33).
Качество получаемых случайных чисел при использовании рекомендуемых пар лагов зависит от величины a: чем оно больше, тем выше размерность пространства, в котором сохраняется равномерность случайных векторов, образованных из полученных случайных чисел. Но в то же время с увеличением лага a возрастает объем памяти, используемой программой, которая реализует датчик. Для практического использования можно использовать такие рекомендации.
Значения (a, b) = (17, 5) хороши для большинства типовых приложений, не использующих векторы высокой размерности со случайными компонентами [9].
Значения (a, b) = (20, 5) применяются в математической системе MATLAB [2], используемой для различных целей, включая ее расширение Simulink для ситуационного моделирования. Автором первой версии этой системы был Д. Каханер [10].
Значения (a, b) = (55, 24) обеспечивают получение чисел, пригодных для применения в большинстве алгоритмов, которые особенно требовательны к качеству случайных последовательностей. В частности, в основном генераторе системы имитацион-
№ 5 (35) 2011
Ü
сС ^
S
и ^
Ü Й
I
00
1 it
со
Si
и
I
s
1 I
i
i
£ §
ного моделирования Actor Pilgrim [6] используются именно эти значения. Одновременно в этой системе применяется и конгруэнтный датчик Р. Шеннона; выбор датчиков осуществляется автоматически.
Значения (a, b) = (97, 33) позволяют получать очень качественные случайные числа и используются в алгоритмах, работающих со случайными векторами высокой размерности [9]. В системе Actor Pilgrim предусмотрена возможность подключения и такого генератора.
Период повторяемости псевдослучайной последовательности фибоначчиевого генератора для чисел с плавающей точкой типа double можно оценить по формуле [10]:
T = (2maxla bl - 1)-2ш, (5)
где ю — число бит в мантиссе вещественного числа.
Следует отметить, что при реализации фибоначчиевого генератора на основе целых 32-разрядных чисел период T получается значительно меньше. В таблице 1 приведены расчетные значения периода генератора в соответствии с формулой (5) для рекомендуемых выше значений лагов при использовании чисел double (ю = 52, 252 = 4503599627370496). Для сравнения в первой строке табл. 1 приведена величина периода генератора Р. Шеннона.
Реализация алгоритма лаг-генератора. Процедура получения очередного x¡ с помо-
щью фибоначчиева генератора требует знать max {a, b} предыдущих сгенерированных случайных чисел. При программной реализации для хранения сгенерированных случайных чисел используется конечная циклическая очередь на базе массива типа double.
Однако для запуска фибоначчиевого генератора сначала необходимо выполнить процедуру начальной загрузки: нужно откуда-то загрузить первые max {a, b} случайных чисел в циклическую очередь, до того как генератор начнет работать. Эти начальные числа могут быть получены с помощью одного из известных конгруэнтных генераторов [8, 15], но лучше всего для этого использовать рассмотренный выше генератор Шеннона.
В качестве примера рассмотрим работу генератора со значениями лагов (a, b) = (17, 5). Первое число, сгенерированное по формуле (4), имеет номер , = 18 при условии, что первые 17 чисел (i = 1, 2, ..., 17) получены конгруэнтным генератором. Циклическая очередь в этом случае имеет следующий вид (рис. 2).
Первоезначение: /'=18 (получено фибоначчиевым датчиком)
b = 5
а = 17
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Первые 17 стартовых значений, сгенерированных конгруэнтным датчиком
к-и
Рис. 2. Генерация первого числа
Таблица 1
Расчетные значения периода программного генератора
Лаги Значение 2' Период t
a b
Шеннона — 5,368709120000000 • 108
17 5 131072 5,902913067590782 • 1020
20 5 1048576 4,722361979270017 • 1021
55 24 4503599627370496 1,622592768292133 • 1032
97 33 158456325028528675187087900672 7,136238463529799 • 1044
104
№ 5 (35) 2011
Намеренно пропустим рассмотрение процедур генерации чисел со второго по пятое. После генерации шестого числа (( = 23) состав циклической очереди выглядит, как это представлено на рис. 3.
Шесть значений: i = 23 (получены фибоначчиевым датчиком)
Г"
]xi-a % Xi-2 */-3 Xi-4
О
b = 5
1 2 3 4 5 6 7
10 11 12 13 14 15 16 17 а = 17
Остаток стартовых значений,
полученных конгруэнтным датчиком <-►
Рис. 3. Получение шестого числа
После генерации седьмого числа (( = 24) циклическая очередь примет вид, отраженный на рис. 4.
Семь значений: /' = 24 (получены фибоначчиевым датчиком)
1
\xi-a xi -*М Xi-2Ixi-3 Xi-4IXi-5I
12 13 14 15 16 17 b = 5
массива имеет номер 0, а последний — номер а (в нашем случае 17).
Первые х1, х2, ..., ха значений предварительно получены конгруэнтным методом и помещены в элементы массива с номерами 1, 2, ..., а. Этими значениями пользоваться не следует.
Обозначим через к номер элемента массива, в котором появляется сгенерированное число при очередном обращении к фи-боначчиевому генератору. После первого обращения к генератору получается первое число. Его порядковый номер i с учетом предварительно полученных а значений из конгруэнтного генератора всегда равен а + 1 (в нашем примере i = 18). Помещается это число в элемент массива к = 0, откуда оно может быть передано для использования.
Итак, номер ячейки, где получается первое х, к = 0. Последующие значения при очередном обращении к фибоначчиевому генератору получаются по рекуррентной формуле:
k =
1 23456789 10 11 а =17
Остаток стартовых значений,
полученных конгруэнтным датчиком <-►
Рис. 4. Генерация седьмого числа
Из рисунков 2, 3 и 4 видно, как в процессе работы фибоначчиевого генератора начальные числа в циклической очереди, предварительно полученные конгруэнтным методом, постепенно заменятся новыми значениями. С помощью этих же рисунков можно сформулировать следующие правила размещения сгенерированных чисел в циклической очереди, ориентируясь на рекомендуемые значения лагов.
Размещение сгенерированных чисел в циклической очереди. Циклическая очередь — это массив для размещения чисел типа double. Длина массива равна a + 1. В нашем примере a = 17. Первый элемент
k-1, если k > 0 a, если k = 0
(6)
Обозначим через т номер элемента массива, в котором находится текущее значение лага а при очередном обращении к фи-боначчиевому генератору. При первом обращении к генератору это значение лага находится в элементе т = а. Номер ячейки, где находится число х, _ а при первом обращении к генератору, т = а (в нашем примере т = 17). Последующие значения т при обращении к этому генератору получаются по формуле:
m =
m -1, если m > 0 a, если m = 0
(7)
Номер элемента массива, в котором находится текущее значение лага а при очередном обращении к генератору, обозначим через п. При первом обращении к генератору это значение лага находится в элементе п = Ь. Номер ячейки, где находится х_ь при этом первом обращении, п = Ь (в данном
2 §
105
№ 5 (35) 2011
¿а ¡С
S
и ^
i
is
Si
со
и
со
Si
и
ii S
1
t
i
is §
случае п = 5). Далее значения п при работе фибоначчиевого генератора получаются по аналогичной формуле:
Гп -1, если п > 0
п = \ . (8)
| а, если п = 0
Далее с помощью фибоначчиевого датчика по формуле (1) можно получать последовательность псевдослучайных чисел с очень большим периодом Т и очень хорошими статистическими свойствами, используя соотношения, которые показывают изменения местоположений чисел в циклической очереди: х1 — формула (6), хн а — формула (7) и х-ь—формула (8).
Нормальный закон распределения
Обратимся теперь к рассмотрению методов получения последовательностей псевдослучайных величин, используемых при исследовании случайных процессов, связанных с рисками, под которыми будем понимать нежелательными события или явления, следствиями которых может быть ущерб. Часто в качестве базового распределения в таких задачах используется нормальное или гауссово распределение (рис. 5).
и имеет дисперсию D (х) = о2. Соответственно, среднеквадратичное отклонение рав но а. Плотность распределения описывается функцией f (х):
1 . Г (х-ц)2
f (x) =
^exp
(9)
ц = 2
Рис. 5. Плотность распределения нормального закона
Это, несомненно, одно из наиболее важных и часто используемых непрерывных распределений. Оно симметрично относительно математического ожидания М( х) = ц
оТ2Л "г| 2о2
Нормальный закон при решении задач, связанных с рисками, часто считается идеализированным законом распределения. «Хвосты» распределения — крайние положительные и отрицательные ненулевые значения х уходят, соответственно, в плюс и минус бесконечность. При моделировании периодов времени с использованием такого распределения можно получить, в том числе, и отрицательные значения времени. Но если мы в модели хотим анализировать моменты возникновения рисковых событий для прогнозирования, то понятие отрицательного времени может стать абсурдным по следующей причине.
Получение отрицательного значения имеет очень маленькую вероятность, но и рисковые события также маловероятны, поскольку происходят редко. Абсурдность заключается в том, что при какой-то случайной реализации моделируемого процесса теоретически возможен прогноз на прошлое. Конечно, вероятность такого прогноза может быть очень мала. Появление такого прогноза можно считать методической ошибкой и не учитывать — отбраковывать. Но отбраковки приведут к «усечению» функции плотности распределения слева от нуля и к некоторому увеличению матожи-дания по сравнению с ц, т. е. к искаженным результатам, и распределение уже не будет нормальным.
Однако несмотря на «идеальность» нормальный закон является базовым, без которого часто бывает невозможно получить даже неточное решение.
Иногда для получения нормально распределенных величин с плотностью распределения (9) используются генераторы или таблицы, которые выдают симметричное относительно нуля распределение вида
106
№ 5 (35) 2011
f (z) =
ЖexpI-2z
(10)
являющееся нормальным, с матожиданием ц = 0 и среднеквадратичным отклонением a = 1. Для получения случайной величины x с плотностью распределения (9) и с произвольными значениями ц и a через величину z, используется преобразование:
x = ц + a • z.
Поэтому для математического ожидания ц и среднеквадратичного отклонения a используют и другие названия: ц — смещение; a — масштаб.
Применение центральной предельной теоремы. Самый простой и быстрый в смысле времени получения очередного числа генератор основан на центральной предельной теореме. Берется сумма 12 равномерно распределенных на отрезке [0,1] случайных величин, результат уменьшается на 6 (сдвигается влево), и случайная величина, имеющая плотность распределения (10), получена. Генератор работает очень быстро. Вместе с тем, в смысле качества получаемых чисел он имеет следующее нехорошее свойство: «хвосты» плотности распределения сравнительно короткие; они уходят в ноль на расстоянии 6a в положительную и отрицательную стороны от матожидания.
Но такое ограничительное свойство не мешает получать точные результаты в подавляющем большинстве практических оценочных задач, связанных с физикой, химией, техникой, военным делом или управлением: во многих учебниках по теории вероятностей показывается, что разброс даже в 3a является наиболее подходящим. Однако при моделировании рисковых ситуаций такой генератор использовать нельзя, поскольку вероятность того, что случайная величина x находится довольно далеко за пределами интервала [ц - 6a, ц + 6a], не равна нулю.
В системе Actor Pilgrim генератор, работающий на основе центральной предельной теоремы, реализован с помощью программной функции central.
Применение преобразования Бокса- | Мюллера. При моделировании «хвостов» г§ распределения могут возникнуть вопросы: ,2
«А что такое вероятность на "хвосте" рас- ^
1 /
пределения, равная, например, '/2000 при исследовании рисков? Насколько она значима?»
Для ответа рассмотрим пример: эксперт-но оценена вероятность того, что на протяжении всех дней в течение ближайших 5,5 лет в какой-то экономической системе (или на каком-то объекте) не произойдет ни одного рискового события конкретной заданной природы, равна р = 0,9995. Тогда величина
q = 1 - р = 1 -0,9995 = 0,0005 = '/2000.
Есть вероятность того, что хотя бы одно рисковое событие (которых может быть и больше!) той же природы произойдет в один из дней в течение ближайших 5,5 лет, т. е. в ближайшие 2000 календарных дней. И если других данных нет, то можно предположить появление одной из следующих двух ситуаций.
1. Допустим, что жизненный цикл системы находится на этапе установившегося режима работы. В этом случае вероятность qд того, что рисковое событие произойдет в ближайший или любой другой день, — малая величина ^ >> qд), и с позиций теории полезности для риск-менеджера q будет существенно более значима, чем qд. Причем для исследования такой ситуации можно проводить вероятностные оценки типа схемы независимых испытаний Я. Бернулли.
2. Однако если жизненный цикл системы в течение этих 5,5 лет по каким-то причинам выйдет из установившегося режима на реорганизацию или на завершающий этап существования, то вероятность qд может резко увеличиваться на протяжении последовательности дней неизвестной длины. В этом случае оценка разницы величин q и qд становится безразличной, хотя понятно, что q > qд. Опытный риск-менеджер, знакомый с теорией полезности, будет считать обе вероятности q и qд одинаково значимыми на протяжении всех 2000 дней. В разряд
107
№ 5 (35) 2011
ja ¡E
S
и ^
Ü Й
I
со 1
со
Si i s
s 1
i
i
£ §
незначимых обе вероятности можно будет отнести только после удачного завершения всего периода 5,5 лет.
Отсюда следует вывод, что при рассмотрении периода времени длиной в 5,5 лет вероятность q = 0,0005 = 1/2000 достаточно велика для прогнозирования как экономических рисков, так и рисков промышленных аварий или природных катастроф.
Рассмотрим генератор нормального распределения, который основан на применении преобразования Бокса-Мюллера [11]. Он программируется сложнее, чем генератор, основанный на центральной предельной теореме, но обеспечивает высокую точность.
Одна из модификаций этого метода заключается в следующем. Введем в рассмотрение три служебных переменных: v2 и е. Далее используем генератор равномерно распределенных на интервале [0, 1] чисел и получим два случайных числа г1 и г2. После этого вычисляются три величины:
^ = -1 + 2г1, у2 = -1 + 2г2 и в = + у22.
Если в > 1 или в = 0, то получается новая пара чисел г1 и г2, и цикл начинается снова. Но если выполняется неравенство 0 < в < 1, то вычисляются два числа и z2, распределенные по нормальному закону, вид которого определяется формулой (10):
-2ln s
S
—2ln s
s
Далее в отношении чисел z1 и z2 применяется преобразование
х1 = ц + а • z1 и х2 = ц + а • z2,
и два новых нормально распределенных числа х1 и х2 можно считать полученными.
Рассмотренная модификация преобразования Бокса-Мюллера имеет следующее эксплуатационные характеристики:
1) необходимы 127 пар случайных чисел г1 и г2 для получения 100 пар нормально распределенных случайных чисел х1 и х2 [15];
2) «хвосты» плотности распределения, получаемые с помощью генератора,
уходят на «невероятно» большое расстояние от матожидания ц, практически в плюс и минус бесконечность, приблизительно на 1,7 х 10308 х а — максимальное значение числа типа double, помноженное на величину а.
В системе Actor Pilgrim этот генератор реализован с помощью функции gauss.
Асимметричные законы распределения
Важнейшим показателем, характеризующим случайную величину, является ее математическое ожидание (или среднее значение). Реже используются такие характеристики как:
• мода — наиболее вероятное значение;
• медиана — значение, которое случайная величина превышает с вероятностью 50%.
Мода и медиана на практике редко используются отчасти потому, что их значения иногда совпадают с математическим ожиданием — особенно в технике и машиностроительных технологиях. Однако в экономике исключения из этого «правила» встречаются довольно часто. Дело в том, что величины, которые принимают только неотрицательные значения (цена, время и т. п.), не могут подчиняться нормальному закону распределения. Такие случайные переменные, как правило, подчиняются законам распределения вероятностей с ярко выраженной асимметричностью, и их математические ожидания, моды и медианы не совпадают.
Типичными примерами таких законов являются логнормальное, экспоненциальное, эрланговское, гамма- и гиперэкспоненциальное распределения. У них есть два общих свойства:
• распределения справедливы только для неотрицательных значений x;
• функции плотностей распределений имеют правые «бесконечные хвосты», асимптотически сходящиеся к нулю, что дает возможность моделировать редкие маловероятные события, происходящие при большом превышении случайной величиной x значений матожидания, моды и медианы, а это
108
№ 5 (35) 2011
Таблица 2 §
Примеры ассиметричных распределений случайных величин
Случайные величины исследуемого экономического процесса Типовое распределение
Стоимость акций при невозможности банкротства [14] Логнормальное
Изменение цены базового актива в моделях ценообразования (свойство волатильности) [12] Логнормальное
Время ожидания страхового события Экспоненциальное
Интервал между приходами заявок в простейшем потоке заявок на обслуживание Экспоненциальное
Время продажи (при наличии спроса) одной единицы товара Экспоненциальное
к единиц товара2 Эрланга, к > 1
Интервал между приходами заявок в групповом потоке Гамма, 0 < к < 1
Интервал между заявками в идеализированном групповом потоке Гиперэкспоненциальное, 0 < к < 1
важно при исследовании процессов, связанных с рисками.
В таблице 2 дано несколько примеров случайных величин, которые обычно распределены асимметрично.
Логнормальный закон распределения. Рассмотрим учебный пример, числовые данные для которого заимствованы из справочника [14]. Допустим, что коммерческая организация имеет некую денежную сумму и должна ею выгодно распорядиться, для чего может выбрать один из двух вариантов:
1) поместить деньги на депозитный счет конкретного банка, который начисляет на вклады 24% годовых;
2) приобрести акции некого эмитента — надежной компании.
Первый вариант. Рассматривая будущую стоимость капитала при вложении его в банк как случайную величину с нулевым стандартным отклонением, нетрудно заметить, что по прошествии одного года организация получит 1,24 руб. на каждый первоначально вложенный рубль. Эту сумму выдает банк, и она гарантирована, если счи-
тать риск возникновения форс-мажорных обстоятельств (например, дефолта) невозможным событием. Поэтому 1,24 руб. в данном случае есть одновременно и мода, и медиана, и матожидание.
Второй вариант. Предположим, что сегодня рыночная стоимость акций некоторого эмитента составляет 1 руб. за штуку. Будущая же их стоимость по прошествии одного года является случайной величиной с лог-нормальным распределением, математическим ожиданием в 1,3 руб. за штуку и стандартным отклонением в 0,4 руб. за штуку (рис. 6). Этот вариант соответствует первой строке табл. 2.
* М
шт./руб
2 Назначение и особенности параметра k будут объяснены ниже при рассмотрении эрланговского, гамма- и гиперэкспоненциального распределений.
3,0 X
руб.
Рис. 6. Логнормальное распределение будущей стоимости акций
Поскольку закон распределения асимметричный, значения x для моды, медианы и ма-
109
-N ПРИКЛАДНАЯ ИНФОРМАТИКА
№ 5 (35) 2011 ' -
ja ¡E
S
и ^
Ü Й
I
со 1
со
Si i s
s 1
i
i
e §
тематического ожидания для будущей стоимости акций различны и соответственно равны 1,135, 1,243 и 1,3 руб. за штуку (на рис. 6 они отмечены тремя пунктирными линиями).
На данном примере покажем, что в процессе принятия решения выбор меры для анализа вариантов крайне важен. Например, если, анализируя распределение, сравнивать будущую стоимость нашего капитала для обоих вариантов инвестирования по моде (1,135 < 1,24), то более выгодным вариантом следует считать вложение денег на депозитный счет в банк. Если же сравнивать по медиане (1,243 = 1,24), то оба варианта приблизительно равны по выгодности. Но сравнение математических ожиданий (1,3 > 1,24) выводит вариант инвестирования в акции на первое место по выгодности.
Рассмотрим особенности генерации лог-нормально распределенных величин. Генератор реализуется просто:
1) в качестве базового генератора используется программный генератор, точно реализующий нормально распределенную последовательность чисел, реализованный, например, на основе преобразования Бокса-Мюллера: генерируется нормально распределенное число м с математическим ожиданием ц и среднеквадратичным отклонением а;
2) выполняется возведение в степень х = ем, где е = 2,718281828459045235... — стандартная математическая константа (иррациональное число); результат х берется в качестве логнормальной величины, распределенной с параметрами ц и а.
Обратное преобразование также справедливо: если х — логнормально распределенная случайная величина с параметрами ц и а, то м = 1пх есть нормально распределенная величина с матожиданием ц и среднеквадратичным отклонением а.
Функция плотности логнормального распределения выглядит так:
f (X) =
1
XG
42П
exp
(ln x -ц)2 2о2
X > 0 X < 0
Следует отметить, что ц и а никоим образом не являются матожиданием и среднеквадратичным отклонением логнормально распределенной величины х, а их размерности не соответствуют реальности. На рисунке 6 третья (правая) пунктирная линия показывает матожидание, но это отнюдь не величина ц.
Соответствующие логнормальному распределению математическое ожидание цлн и дисперсия а2лн связаны с базовым нормальным распределением соотношениями:
цлн = ец+0 /2; агпн = (ea - 1)е2ц+° .
Как правило, параметры ц и а заранее неизвестны. Чаще известны реальные цлн и алн — на основе статистики, прогнозов или экспертных оценок. Поэтому, чтобы воспользоваться базовым нормальным распределением, ц и а можно вычислить с помощью следующих уравнений:
ц = 1п(Цлн)-о2/2; о= ln
I_ лн_
1 ехр[21п(Цлн)]
+ 1
В системе Actor Pilgrim есть два программных генератора, реализующих лог-нормальное распределение:
muller — используется, когда известны
Цлн и СТлн;
lognorm — используется, когда известны ц и а.
Итак, рассмотрено одно из важнейших распределений — логнормальное, функция плотности которого имеет асимметричный вид. Однако существуют и другие не менее важные распределения, базовым законом распределения для которых является экспоненциальный закон.
Экспоненциальное распределение занимает очень важное место при проведении системного анализа экономической деятельности. Этому закону распределения подчиняются многие явления, например:
• время ожидания страхового события (третья строка табл. 2);
2
№ 5 (35) 2011
• время поступления заказа на предприятие;
• время между посещениями покупателями магазина-супермаркета;
• продолжительность телефонных разговоров;
• срок службы деталей и узлов в компьютере.
Рассмотрим это распределение подробнее. Если вероятность наступления события на малом интервале времени At очень мала и не зависит от наступления других событий, то интервалы времени между последовательностями событий распределяются по экспоненциальному закону с плотностью вероятностей
1 е^/ц, x > 0 Р(x) = <¡ ц .
0, x < 0
Особенность этого распределения — его математическое ожидание M (t) = ц и дисперсия и D (t) = ц2. Часто вместо матожида-ния используется параметр X = 1/ц, который называется интенсивностью процесса — поступления заявок, обслуживания или какого-то другого. Математическое ожидание экспоненциально распределенных величин равно среднеквадратичному отклонению, что является одним из основных свойств этого распределения:
= ц = 1/X.
Получение экспоненциально распределенных случайных величин выполняется с использованием известного метода обратных функций [2]. В соответствии с этим методом, если сгенерировать случайную величину w, равномерно распределенную на интервале [0, 1], то ее преобразование
-ц-lnw
представляет собой экспоненциально распределенную величину.
В системе Actor Pilgrim этот генератор реализован в виде функции expont.
Распределение Эрланга. Обычно распределение А. Эрланга используется в слу-
чаях, когда длительность какого-то процесса можно представить как сумму к элементарных последовательных составляющих, распределенных по экспоненциальному закону. Если обозначить математическое ожидание длительности каждой элементарной составляющей как ц, то плотность вероятностей распределения Эрланга представляется следующей формулой:
f ( x ) =
x
- x / ц
цк -(k-1)! 0 ,
x>0. x<0
Математическое ожидание случайной величины, распределенной по закону Эрланга, М( х) = к • ц. Дисперсия такого распределения D( х) = к -ц2.
Очевидно, что при к = 1 получаем экспоненциальное распределение. Разновидности этого распределения для разных к > 1 представлены на рис. 7.
к = 1 — экспоненциальное распределение ц = 2
Рис. 7. Распределение Эрланга
Генерация числовых последовательностей, распределенных по закону Эрланга, производится с учетом того, что этот закон совпадает с экспоненциальным при к = 1. Если щ — случайная величина, равномерно распределенная на интервале [0, 1], то величина (-ц- 1п щ) представляет собой экспоненциально распределенную величину. Поэтому случайная величина, распределенная
«
0
1
111
№ 5 (35) 2011
по закону Эрланга с параметрами k и ц, определяется выражением
(- ln
W:
где w¡ — независимые случайные переменные, равномерно распределенные на интервале [0, 1]; i = 1, 2, ..., k.
В системе Actor Pilgrim этот генератор реализован в виде функции erlang.
Гамма-распределение. Появление групповых потоков в сложных экономических системах приводит к резкому увеличению средних длительностей различных задержек (заказов в очередях, задержек платежей и др.), а также к увеличению вероятностей рисковых событий или страховых случаев. Для моделирования таких потоков используются гамма- и гиперэкспоненциальное распределения.
Функция плотности гамма-распределения определяется следующей формулой:
ja ¡E
S
и ^
I Й
I
со 1
со
Si
и
S
S 1
i
i IS
с §
f (x) =
x
- x / Ц
-r(k) 0 ,
x > 0 x < 0
3 В математической литературе часто вместо латинской буквы G при обозначении гамма-распределения пишут греческую Г, но в данном случае буква Г занята для гамма-функции Л. Эйлера (Прим. автора).
112
то Г(к) = (к -1)!, и получается распределение А. Эрланга. Соответственно, при к = 1 получается экспоненциальное распределение.
Перейдем к рассмотрению генерации гамма-величин. Для этого потребуется знать два следующих свойства гамма-распределения.
Свойство 1. Параметр ц называется коэффициентом масштабирования, поскольку гамма-распределение обладает следующим свойством: если имеют место распределения х ~ й(к,ц) и у ~ б(к,1), то справедливо масштабное соотношение х = ц- у.
Свойство 2. Параметр к в общем случае состоит из суммы целой и дробной частей: к = с + d, где с — целая часть; d — дробная часть. Поэтому необходимо будет учесть такое свойство: если имеют место распределения х ~ б(к1,ц) и у ~ й(к2,ц), то
(x + у)
- G(k2 + k2
|)
где Г(к) — гамма-функция Л. Эйлера.
Гамма-распределение случайной величины х с параметрами к и ц обозначается3 как х ~ й(к, ц). Параметры к и ц — вещественные величины, к > 0, ц > 0. Математическое ожидание и дисперсия случайной величины х определяются соотношениями М (х) = к • ц и D (х) = к • ц2.
При моделировании систем массового обслуживания групповой поток заявок получается тогда, когда интервал времени между поступающими на обслуживание заявками имеет гамма-распределение с параметром к < 1. Если параметр к — целое число,
В нашем случае к1 = с и к2 = d.
Учитывая свойство возможности масштабирования по параметру ц, указанное выше, достаточно получить гамма-величину для ц = 1. Переход к другим значениям параметра ц осуществляется простым умножением. Генерация каждой гамма-распределенной величины х выполняется в три этапа.
Этап 1. Вычисляется гамма-распределен-ная величина для с — целой составляющей. Используется тот факт, что распределение й (1,1) совпадает с экспоненциальным. Если щ — случайная величина, равномерно распределенная на интервале [0, 1], то (- 1п щ) ~ 6(1,1). Гамма-величина целой части распределена по закону Эрланга с параметрами с и 1. Поэтому
Х (- ln
w
G(c,1),
где щ — независимые случайные переменные, равномерно распределенные на интервале [0, 1]; i = 1, 2, ..., с.
Этап 2. Осталось вычислить гамма-рас-пределенную величину для d — дробной составляющей к. Формул для получения d
№ 5 (35) 2011
не существует, но есть несколько итерационных способов точного расчета этой величины. Рассмотрим алгоритм без доказательства (рис. 8). В нем использованы следующие переменные и константы: а — гамма-величина дробной составляющей к: а ~ G ^,1); d — дробная часть к;
у и у — пара независимых случайных величин, распределенных на [0,1]; Ь и у — рабочие переменные; e = 2,718281828459045235 — стандартная математическая константа.
В результате выполнения алгоритма получена гамма-величина а дробной составляющей к: а ~ G (d ,1).
Этап 3. После сложения целой и дробной составляющих и перемножения суммы на масштабный множитель ц получаем значение новой гамма-распределенной величины х:
a + £ (- ln w
i=1
G(k, ц).
В системе Actor Pilgrim соответствующий программный генератор реализован с помощью программы gamma.
Гиперэкспоненциальное распределение. Рассмотренный выше генератор гамма-распределения обеспечивает достаточную точность, но имеет эксплуатационный недостаток: он работает не очень быстро из-за сложных итерационных вычислений при нецелых значениях k, в том числе при k < 0 для моделирования групповых потоков. Однако исследователи, знакомые с теорией массового обслуживания, знают, что задержка заявок в очереди, если входной поток не пуассоновский (в том числе групповой), в основном определяется первыми двумя моментами случайного интервала поступления заявок в очередь, а не видом закона распределения. Поэтому, если необходимо только получить групповой поток при известных математическом ожидании и дисперсии этого интервала для моделирования очереди, а дополнитель-
Í Выделить &.
I к = с + d
со
0
1
а = [vj v)
v. - V
а -1 In 1
1-у
1
b = 2 a
Рис. 8. Алгоритм вычисления гамма-величины дробной составляющей
ные тонкости распределения интервала поступления и формирования группы заявок не входят в интересы исследователя, то рекомендуется использовать не гамма-, а гиперэкспоненциальное распределение. Соответствующий генератор работает намного быстрее и несильно влияет на время выполнения основной моделирующей программы.
113
№ 5 (35) 2011
Плотность гиперэкспоненциального закона распределения определяется формулой
п
1(х) = Х Pi(х) • 4хр/(У ),
/=1
где п — общее число экспоненциально распределенных величин; 1ехр / ( У ) — функция распределения случайной величины у по экспоненциальному закону;
р,(х) — вероятность того, что в составе распределения 1(х) «проявляется» распределение 4хр / ( У/); / = 1, 2, ..., п.
Соответственно справедливо выражение:
—е-у1 у: > 0 Ц/ 1 .
f (У) =
exp/\У / >
о
ja ¡E
S
u ^
i й
I
со
I it
00 Si
H
s
s 1
i
i
e §
у I < о
Закон распределения назван гиперэкспоненциальным в связи с тем, что математическое ожидание случайной величины х меньше величины среднеквадратичного отклонения. Первые два момента этого закона определяются по формулам
п
М( х) = ^цр, — первый момент — мате-
=1
матическое ожидание;
п
М( х2) = 2^ц2 РI — второй момент.
=1
Генератор для формирования группового потока можно получить, используя только два слагаемых: п = 2. Обозначим через константу к вероятность р1 (х), а через 1 - к обозначим р2 (х), 0 < к <1. Тогда формула плотности распределения для работы этого генератора имеет вид:
f (x) = k—ev * + (1 - k) ■ где yi > 0, У2 > 0.
lim
^2
- У2 / ^2
Обозначим ц = ц и получим два момента и дисперсию:
М( х) = к • ц — математическое ожидание; М( х2) = 2к • ц2 — второй момент; D(х) = (2к - к2) •ц2 — дисперсия.
Сопоставление формул для M (x) и D (x) показывает, что при k < 1 среднеквадратичное отклонение превосходит математическое ожидание.
В результате применения такого датчика генератор заявок создаст поток со следующими свойствами: по истечении интервала времени, распределенного по экспоненциальному закону со средним значением ц, из генератора однократно выйдет группа заявок, между которыми нулевой интервал. Случайный размер группы фактически определяется схемой независимых испытаний Я. Бер-нулли. Математическое ожидание n числа заявок в группе вычисляется по формуле
1
n = —, k
откуда следует, что средний размер группы обратно пропорционален параметру распределения k.
Пример 1. При k = 0,5 получим n = 2. Отношение дисперсии к математическому ожиданию в квадрате D( x)/ M( x)2 = 3, т. е. среднеквадратичное отклонение превышает ц в V3 ~ 1,73 раза.
Пример 2. При k = 1 групповой поток невозможен, поскольку получим экспоненциальный закон распределения, группа будет состоять из одной заявки, а D( x) / M( x)2 = 1.
Программный генератор действует следующим образом. При его работе возникает необходимость генерации двух взаимно независимых величин w1 и w2, распределенных равномерно на интервале [0,1]. При каждом обращении генератор выполняет такие действия:
1) генерируется псевдослучайная величина w1;
2) если w1 < k, то генерируется величина w2, и в качестве очередного значения x выдается -ц- In w2;
3) иначе, при w1 > k в качестве x выдается значение 0.
Понятно, что время выполнения таких действий очень мало'.
В системе Actor Pilgrim соответствующий программный генератор для генерации рас-
114
№ 5 (35) 2011
смотренного группового потока обеспечивается программой hyperexp.
Реализация законов распределений в системе имитационного моделирования Actor Pilgrim
Характерным примером системы, предназначенной для моделирования и исследования экономических процессов, является Actor Pilgrim [5, 6]. Она получила распространение в экономических вузах в последнее десятилетие. В ней для каждого процесса, входящего в состав имитационной модели, обеспечивается независимая последовательность псевдослучайных чисел, т. е. в рамках одной имитационной модели фактически создаются виртуальные независимые генераторы.
Все законы распределения случайных величин, включая любые новые распределения, вводимые пользователем — разработчиком моделей, получаются посредством двух базовых программных генераторов, обеспечивающих равномерное распределение:
• фибоначчиевого лаг-генератора с лагами (55, 24);
• конгруэнтного генератора Р. Шеннона, который используется как для начальной загрузки лаг-генератора, так и в каче-
стве самостоятельного быстродействующе- | го генератора с хорошими статистическим г§ характеристиками, если к псевдослучайной ¿g последовательности не предъявляются тре- ^ бования по обеспечению супердлинного пе- ^ риода.
В фибоначчиевом генераторе пара лагов (55, 24) применяется по умолчанию. Переключение на другую пару лагов (97, 33) производится изменением настроечных параметров системы. Практика показывает, что необходимости в таком переключении не возникает из-за очень хороших свойств генератора с лагами (55, 24).
При проведении имитационного моделирования повышенные требования к генератору предъявляются в модельных процессах типа:
actor — генератор акторов с бесконечной емкостью;
produce — управляемый процесс (непрерывный, пространственный);
serve — узел массового обслуживания с произвольным числом параллельных каналов, с возможностями прерываний.
С помощью комбинаций именно этих процессов в имитационных моделях выполняются следующие действия:
• моделируются редкие рисковые события на длинной оси времени, при одновременном моделировании хозяйственно-
Таблица 3
Распределения случайных величин
Параметр Программа Реализуемый закон распределения
cent central Нормальный закон: центральная предельная теорема
norm gauss Нормальный закон: точная реализация
logn lognorm Логнормальный закон: задается параметрами ц и о
muln muller Логнормальный закон: задается через цлн и олн
expo expont Экспоненциальное распределение
erln erlang Распределение Эрланга
gamm gamma Гамма-распределение
hype hyperexp Гиперэкспоненциальное: для группового потока
№ 5 (35) 2011
экономической, финансовой, военной или иной деятельности, для которой эти события представляют угрозу;
• при моделировании проявляются «хвосты» распределений вероятностей рисковых событий практически бесконечной длины, что дает возможность для составления прогнозов.
Поэтому для моделирования таких процессов используются только фибоначчие-вые датчики. Процессы других типов используют конгруэнтный генератор. Выбор генератора осуществляется автоматически на основе анализа типа процесса, из которого пришло обращение для получения случайной величины. Даже если пользователь — разработчик моделей захочет в явном виде обратиться к базовой программе получения псевдослучайных величин, он будет использовать обращение к функции random, которая сама вызовет либо программу kahaner (фибоначчиев лаг-генератор), либо программу shannon (конгруэнтный генератор). Все перечисленные программы реентерабельные чтобы обеспечить каждый == моделируемый процесс своей независимой £ последовательностью чисел. 4с Система Actor Pilgrim обеспечивает ! стандартные программы для автоматиче-<3 ского получения рассмотренных выше за-^ конов распределения случайных величин l (табл. 3), которые подключаются автоматиЦ чески (следует отметить, что в Actor Pilgrim 5 имеются стандартные программы для полу-| чения и других распределений). I Стандартный способ получения случай-si ных величин, распределенных по нужному I закону, состоит в том, чтобы в моделирую-I щих функциях actor, produce или serve под-ss ставить соответствующий параметр, указан-^ ный в первом столбце табл. 3. Иногда для § каких-то целей пользователю-модельеру Е| потребуется из произвольного места ими-|| тационной модели, из любого ее узла не-| зависимо от его типа обратиться к конкрет-sg ному генератору-программе. В этом случае он сможет в явном виде использовать со-¿s ответствующую программу, задав при об-
ращении к ней необходимые аргументы [5]; имена программных генераторов перечислены во втором столбце табл. 3. Кроме того, предоставляется возможность подключения пользовательских генераторов, реализующих любое распределение, если эти генераторы используют базовую функцию random.
Заключение
Выше рассмотрены теоретические основы и способы реализации программных генераторов, использующих алгоритм Фибоначчи с лагами запаздывания и конгруэнтный метод, обеспечивающих высокие статистические характеристики следующих генерируемых псевдослучайных последовательностей:
• двух базовых генераторов, обеспечивающих получение равномерно распределенных величин на интервале [0,1], основного лаг-генератора kahaner и конгруэнтного shannon;
• семи генераторов, особенностью которых являются конкретные требования, предъявляемые к получению длинных, практически бесконечных «хвостов» распределений в задачах, связанных с анализом рисковых процессов: gauss, lognorm, muller, expont, erlang, gamma и hyperexp;
• быстродействующего генератора нормального распределения central для имитационного моделирования производственных процессов в промышленности, обеспечивающего конечные «хвосты» с длиной 6о.
Корректность используемых генераторов подтверждается как их тестированием, так и многолетней эксплуатацией имитационных моделей экономического назначения.
В Московском финансово-промышленном университете «Синергия» пакет Actor Pilgrim, оснащенный такими программными генераторами, успешно используется в учебном процессе при подготовке специалистов с высшим образованием по направлениям:
«Прикладная информатика» — бакалавры, магистры и информатики-экономисты;
№ 5 (35) 2011
«Математическое обеспечение и администрирование информационных систем» — бакалавры, магистры и математики-программисты.
В учебных планах предусмотрены учебные дисциплины по компьютерному моделированию и основам теории рисков — в сумме 128 часов: 32 часа — лекции, 32 часа — семинары и практические занятия, 64 часа — самостоятельная работа, в процессе которой выполняется курсовой проект в среде Actor Pilgrim.
При выполнении курсовых проектов используется «бригадный» подход, суть которого заключается в том, что несколько студентов изучают конкретный экономический рисковый процесс, его особенности, проводят его структурный анализ и создают послойную имитационную модель этого процесса.
После верификации и калибровки модели с ее помощью проводится экстремальный компьютерный эксперимент, в процессе которого получаются неизвестные распределения рисковых событий с длинными «хвостами», причем далеко не все они совпадают с известными в математической литературе законами.
Конечной целью эксперимента является оптимизация параметров, влияющих на стоимость мероприятий по управлению рисками и возможные потери при наступлении рисковых событий. Таким образом, в учебном процессе студенты приобретают не только соответствующие знания по компьютерному моделированию, но и дополнительные компетенции в области экономических рисков и риск-менеджмента.
Список литературы
1. Девятков В. В., Кобелев Н. Б., Емельянов А. А., Половников В. А., Плотников А. М. Имитационное моделирование как основной способ поддержки принятия решений в современном мире. Об организации имитационных исследований в России. — В кн.: Третья Всероссийская научно-практическая конференция по имитационному моделированию и его применению в науке и промышленности «Имитационное моделиро-
вание. Теория и практика». Т. 1. СПб.: Санкт-Пе- §
"5
тербург, ЦНИИ ТС, 2007. С. 37 - 47.
2. Дьяконов В. П. MATLAB 6.5 SPl/7 + Simulink g 5/6. Основы применения. Серия «Библиотека ** профессионала». М.: СОЛОН Пресс, 2005. — ^ 800 с.
3. Емельянов А. А. Имитационное моделирование в управлении рисками. СПб.: Инжэкон, 2000. — 376 с.
4. Емельянов А. А. Имитационное моделирование экономической динамики // Прикладная информатика. 2010. № 1 (25). С. 105 - 118.
5. Емельянов А. А, Власова Е. А, Дума Р. В., Емельянова Н. З. Компьютерная имитация экономических процессов / под ред. А. А. Емельянова. М.: Маркет ДС, 2010. — 464 с.
6. Емельянов А. А., Власова Е. А. Имитационное моделирование экономических процессов. — В кн.: Четвертая Всероссийская научно-практическая конференция по имитационному моделированию и его применению в науке и промышленности «Имитационное моделирование. Теория и практика». Т. 1. СПб.: Санкт-Петербург, ЦТСС, 2009. С. 19 - 26.
7. Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение. М.: Мир, 1998. — 575 с.
8. Кельтон В., ЛоуА. Имитационное моделирование. Классика CS. СПб.: Питер, BHV, 2004. — 848 с.
9. Кнут Д. Э. Искусство программирования. Т. 2. Получисленные алгоритмы. М.: Вильямс, 2007. — 832 с.
10. Метод Фибоначчи с запаздываниями (Lagged Fibonacci Generator). Википедия, 2011. URL: ru.wikipedia.org.
11. Преобразование Бокса-Мюллера. Википедия, 2011. URL: ru.wikipedia.org.
12. Смывин А. Ю. Модели определения стоимости и управления риском опционных контрактов. Канд. дисс. (08.00.13). М.: Московская финансово-промышленная академия, 2010. — 156 с.
13. Стандарт IEEE 754 - 2008. Википедия, 2011. URL: ru.wikipedia.org.
14. Токарев С. С. Справочник экономиста. М.: Корпоративный менеджмент, 2009. URL: http://www.cfin.ru/finanalysis/math/crook.shtml.
15. Шеннон Р. Е. Имитационное моделирование систем: наука и искусство. М.: Мир, 1978. — 420 с.
