Проблемы генерации случайных величин с фрактальными распределениями Текст научной статьи по специальности «Математика»

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 3 (113) 2012

УДК 519.2:004/421.5:004.7

В. Н. ЗАДОРОЖНЫЙ О. И. КУТУЗОВ

Омский государственный технический университет

Санкт-Петербургский государственный электротехнический университет «ЛЭТИ» им. В. И. Ульянова (Ленина)

ПРОБЛЕМЫ ГЕНЕРАЦИИ СЛУЧАЙНЫХ ВЕЛИЧИН С ФРАКТАЛЬНЫМИ РАСПРЕДЕЛЕНИЯМИ

Вскрываются специфические проблемы генерации случайных величин, описываемых распределениями с тяжелыми хвостами. Обосновывается необходимость использования для их генерации многоразрядных ЭВМ.

Ключевые слова: случайные числа, системы с очередями, моделирование.

Работа выполнена при поддержке гранта РФФИ 12-07-00149-а.

Введение. Фрактальные, т.е. степенные или асим-птотически-степенные [1] распределения — это распределения с тяжелыми хвостами (РТХ). В книге [2] в числе перспективных моделей теории массового обслуживания упоминаются «модели с «тяжелыми хвостами» распределений, характеризующих входящий поток и процесс обслуживания».

В настоящей статье рассматриваются специфические проблемы имитационного моделирования (ИМ) систем с очередями, обусловленные необходимостью генерации случайных величин (сл.в.) с РТХ.

Проблема корректной реализации фрактальных распределений. При ИМ реализация фрактальных распределений сопряжена со специфическими трудностями вычислительного характера. Поясним их на примере распределения Парето Ра(К,°а), описываемого плотностью вероятностей (п.в.)

ctKa

,

(1)

где а — параметр формы, К > 0 — минимум случайной величины х (масштабный параметр). Соответствующая функция распределения (ф.р.) имеет вид:

F(t) = 1 — (K/t)“, t > K

(2)

Из (1) к-й начальный момент распределения Парето (РП) находим в виде:

M(x*) = jW(f)df =

(3)

<хК

a >k,

a-k qo , а < k,

тальных сл.в. обнаруживаются при их практическом использовании в ИМ. Например, в [3] описаны статистические эксперименты с генераторами паретов-ской сл.в. хеРа(К, а), и для некоторых а приведены результаты численного анализа свойств генерируемых величин. Эти эксперименты и численный анализ показывают, что конечная разрядность вещественных чисел ЭВМ приводит к необычайно высоким погрешностям реализации моментов сл.в. х. А эти погрешности, в свою очередь, приводят к серьезным ошибкам при ИМ очередей. Поэтому при разработке общих методов ИМ очередей важно исследовать погрешности реализации моментов РТХ.

Найдем погрешности реализации м.о. сл.в. хеРа(К, а) для всех а>1. Будем рассматривать сл.в. х как идеальную (математическую) сл.в., принадлежащую распределению Ра(К, а), и, соответственно, имеющую в точности все те свойства, которые им определяются, в частности — моменты (3). Величина, реализуемая компьютером, отличается от х, и мы обозначим ее через х.

Точная формула генератора сл. в. хеРа(К, а) дается методом обращения функции (2) и имеет вид х = К/*К^г , или, проще,

x = K/*[z ,

(4)

Отсюда при а>1 находим конечное среднее (м.о.) М(х) = аК/(а—1), при а>2 — конечную дисперсию D(x) = M(x2)-M2(x) = аК2/(а -1)2/(а -2).

Погрешности реализации математического ожидания РТХ. Особые свойства датчиков фрак-

где zeR(0, 1) — базовая сл.в., (БСВ) — непрерывная, независимая, равномерно распределенная в промежутке (0, 1) величина. На ЭВМ БСВ z представляется вещественной в компьютерном смысле, но дискретной в действительности псевдослучайной величиной z, которая принимает равновероятно одно из значений, расположенных между 0 и 1 с малым шагом дискретизации е. Так, в MS Excel е=10-15, а в GPSS е=10-6. Соответствующее формуле (4) преобразование БСВ z дает дискретную сл.в.

(если при этом возможно значение z = 0 , то оно не используется).

Положим, ze{e,2e,3e,...,Ne} , где Ne=1, и найдем

м.о. М(х) = М(К/л/5) при условии, что все значения z,=!e = z'N_1 дискретной сл.в. z равновероятны:

N 1 N І(

мй = ^,р &) = ^І-£= =

і=і N і=і с^г.

N .

= Ш-1'£2І “ ^ЛП^ЇЛГ1) “ =

1=1 1=1

N _і і 1_і N _1

=Ш~1,£і аЫа=КЫа “ =

Выполняя в (11) замену р = 1/а и отбрасывая в правой части малое слагаемое R2 и член с М-1, получаем для Ам более простые асимптотические по N границы:

—1

К

№

а-1

—і

■ < А]уі < К'

аЫ°_

а-1

(12)

= КЫ р-^і-р = КЫ р-1Бы (р) /=1

(а>1, р = 1/а, 0<р<1).

(5)

или, компактнее:

Ш м <Ам<МЫ м,

(13)

Оценим сумму из (5) ^(Р) = 1Гр снизу и сверху,

1=1

используя интегралы:

, (N>1).

1 1

Вычисляя интегралы, получаем:

(* + ГР-1 < 5 (р) < ^~Р~Р

Выполняя здесь разложение функции (N+1)^р в знакочередующийся ряд

(6)

(ЛТ + 1)1_р = АГ1"'’ + (1 -р)ЛЛр + (1 Р*ДГрЧ +... (7)

(N+1)^р=№-р+^-р^-р^2, (8)

где R2>0, и, как видно из (7), R2<0,5(1— р)р^ р-1. Таким образом, R2 составляет в (8) менее №й доли слагаемого (1-р)М-р, которое, в свою очередь, составляет менее N й доли слагаемого №-р. Заметим, что нас интересуют значения N>10®. Подставляя (8) в (6), получаем оценку:

. (9)

1-р 1-р

Наконец, умножая (9) на К№-1, получаем, с учетом (5), оценку для М(х):

±М+^-1<мй<х^^. (10)

1-р

1-р

Сравнение верхней оценки для с М(х) =

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

м М(х) 1 -рАГ^1

При всяком N, если р^1 (а^1), то коэффициент ум^“. Из (10) для абсолютной погрешности Дм = |М(х) - М(х)| = К/(1 - р) - М(х) находим:

КЕ^<Ам<К

1-р

мр (1 + Я2)_ш-1. (11) 1-р

где N^1», М = М(х) = Ка/(а-1).

Как видно из (13), погрешность реализации м.о. паретовской сл.в., обусловленная дискретностью «вещественных» чисел в ЭВМ, с ростом N сходится к нулю (что естественно), но сходится она по асимптотически степенному закону (что в общем случае слишком медленно).

Рассмотрим, например, реализацию сл.в. х при К=1, а =1,1 и N=10®. Значение М(х), согласно (3), в этом случае составляет М(х) = аК/(а-1) = 11. Однако его реализация в ЭВМ занижена на величину Ам. Из (13) находим, что эта величина лежит в пределах от 2,848 до 3,133, т.е. составляет около трех.

Точные расчеты, проведенные в следующем разделе статьи, показывают, что вместо м.о. 11 в рассмотренном случае реализуется м.о. 8,0297. При а= = 1,05 или а =1,01 отличия реализованных м.о. от истинных становятся многократными, и даже двойная точность вычислений (случай N=10“) положения не спасает.

Числовые примеры. В качестве примеров в табл. 1 приведены значения погрешностей Ам и коэффициентов занижения ум при некоторых ае(1, 2) и К=1. Значения М(х) рассчитываются по формуле (5) через сумму SJY(p). При N<10® сумма SJY(p) легко вычисляется на персональной ЭВМ непосредственно. При N>>10® непосредственное вычисление этой суммы затруднительно, и она рассчитывается как сумма двух частей:

т N

5л,(р) = ^Гр + 2]г_р , где т=106. Первая часть, 1=1 1=01+1

т.е. сумма

вычисляется непосредственно.

N

Вторая часть — сумма

вычисляется как сред-

нее арифметическое (д + G)/2 ее нижней д и верхней G оценок, где

N+1

д= Р(Ц =

т+\ N

С=\ГР(Н =

(АГ + 1)1~р-(ді + 1)1~р 1-р

М'-р -т'-Р

1-р

(14)

Такой расчет второй части суммы Бы(р), как показывают численные проверки (выполненные для 0,01<р<0,99), обеспечивает при N>10® получение минимум семи точных значащих цифр. Так, при т=10®, N=10“, р = 0,8 получаем д=4920,75532,

N

= 4

С=4920,75534,

= 4 920,755 33.

Из табл. 1 видно, что даже при N=10“ (т.е. при двойной точности вычислений) погрешности реализации паретовской сл.в. для цифровой ЭВМ фено-

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 3 (113) 2012 ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 3 (113) 2012

1,0E+02 1,0E+00

i,oe-o2e

1,0E-04 1,0E-06 1,0E-08 1,0E-104

Рис. 1. Зависимость погрешностей Дм от N при различных ае(1, 2)

а=1,1 N

+00 1E+0^*. 1E+35 1E+42 1E+49

>0=1,25

a=1,5 "V ' іл.лл -0,0909 ^y = 10,429x

л ллпс. -0,3333 >sv y = 2,4476x -0 2 y = 4,4375x0,2

Рис. 2. Зависимости математических ожиданий и дисперсий от а li

Рис. 3. Зависимости м.о. и дисперсий от а во всем диапазоне а > 0. Ряды 1, 2 - соответственно м.о. и дисперсия при N = 106,

Ряды 3, 4 - соответственно м.о. и дисперсия при N = 1015

менально высоки, особенно при значениях а, приближающихся к единице.

Графики на рис. 1 представляют зависимость погрешностей Ам от N.

Полученные численными методами уравнения регрессии на рис. 1 степенные, достоверность аппроксимации R2 максимальна (равна единице). Показатели степени и коэффициенты уравнений регрессии соответствуют аналитическому результату (13).

Погрешности реализации высших начальных моментов. Используя РП в качестве типичного представителя распределений со степенными и асимпто-тически-степенными хвостами, рассмотрим погрешности реализации его моментов порядка к>2. У РП к-й начальный момент ^к=М(хк) при а<к бесконечен (3) и равен аКк/(а-к) при а>к. Для определения реализованного момента =М(х4) по аналогии с формулой (5) получаем:

Mk=M(xk) = fjx^(xi) = ±-fj

К

= KkN~1^zia=KkN-l^4(iN~1) “ =

1=1 *= 1

]у к к к ^ д г к

= KkN~1'£ji~aNa =KkNa~ =

І=1 /=1

N

= KkNp~1^rp = KkNp~1SN(p}

І=1

(а>к, p = k/a , 0<p<1).

(15)

Чтобы найти границы для значений At, достаточно умножить установленное для суммы SN(p) соотношение (9) на KkNp-1. Выполняя это, находим:

Kkl-NP^+R^+KkN-<fik<Kkl-PNP'\ (16) 1 -р 1 -р

Сравнение верхней оценки для At с цк=аКк/ /(a — k) = Kk/(1—k/a)=Kk/(1 — p) показывает, что реализованный k-й момент At всегда меньше, чем /лк. Из (16) по аналогии с (11), (12) вытекают асимптотические границы для абсолютной погрешности Дfit =Khl(\-р)-ftt реализации момента /лк:

Рис. 4. Зависимости коэффициента вариации от а; а>0. Ряд 1 — при N=106 ряд 2 — при N=10^

-АЛТ“_1<Д^<^“~1. (17)

а

Формула (13), найденная для ц1 = М, — это частный случай формулы для цк (17).

Соотношение (17) устанавливает факт асимптотически степенной сходимости к нулю погрешностей А^к с ростом N, и точно определяет показатель степени (к/а-1), характеризующий скорость сходимости.

Погрешности реализации дисперсии. Из (13) и (17) вытекает, что погрешность реализации дисперсии D(x) паретовской сл.в. с ростом N также сходится к нулю лишь с асимптотически степенной скоростью. Действительно, поскольку дисперсия D(x) = M(x2) — M2(x), а реализуется дисперсия D(X) = М(х2) - М2(х), то погрешность реализации дисперсии составляет

Дс = Б(х) - О(х) = М(х2 )-М2(х)-[М(х2 ) - М2 (х)] = = М2~М2~ [М2(х) - М2(х)]=

=Л//2 - [М(х) - М(х)] • [М(х) + М(х)]=

=А/л2 - Ам • [М(х) + М(х)],

(а>2). (18)

В последнем выражении обе величины А^2 и Ам сходятся к нулю со степенной скоростью, а сумма сходится к константе 2М(х). Отсюда следует, что сходимость к нулю погрешности АD является асимптотически степенной.

Кроме того, в (18) в последнем выражении скорость сходимости к нулю для А^2 характеризуется показателем степени (2/а-1), а для Ам — показателем (1/а — 1). Из этого вытекает, что скорость сходимости к нулю погрешности АD характеризуется показателем (2/а-1), который больше и определяет меньшую скорость сходимости, чем показатель (1/а — 1). Из этого вытекает также асимптотика А^А^2 при N^1».

В табл. 2 приведены погрешности АD и коэффициенты занижения дисперсии ув = В(х)/Ъ(х) для некоторых ае (2, 3) при К=1.

Рассчитанные численными методами зависимости математических ожиданий и дисперсий от а представлены на рис. 2.

Сплошные линии на рисунке — это графики действительных моментов, линии с белыми маркерами — графики моментов, реализуемых при N=10®, и линии с черными маркерами — графики моментов, реализуемых при N=10“

Реализация бесконечных м.о. и дисперсий. У паретовской сл.в. х при а<к момент цк = »; моменты реализованной численно сл.в. х при а<к, естественно, конечные, и формула (15) позволяет их вычислить. Для таких вычислений можно воспользоваться и более простыми оценками (16). Однако еще более простые оценки (17) для случая а<к неприменимы, т.к. при их выводе использовалось соотношение (3), справедливое лишь при а>к.

Вместе с тем очевидно, что абсолютная погрешность реализации моментов и коэффициенты их занижения при а<к бесконечны, поскольку при а<к вместо истинных бесконечных моментов реализуются моменты конечные.

Результаты вычисления моментов, реализуемых при а<к, представлены графиками на рис. 3. Эти графики позволяют составить представление об изменении реализованных моментов во всем диапазоне значений а и при различных e = N—1. Кроме того, поскольку при анализе очередей весьма важной характеристикой распределений является коэффициент вариации, то ниже представлены также и графики изменения его реализации С = ^О(х) /М(х) (рис. 4). На основе соотношения (16) аналитическими методами (выводом аналогичного ему соотношения для коэффициента вариации и взятием соответствующего предела) можно показать, что при а^0 коэффициент вариации С сходится к у[ы . И на рис. 4 такую сходимость можно видеть.

Моделирование вычислительных сетей с фрактальным трафиком. Основной особенностью трафика современных телекоммуникационных сетей с коммутацией пакетов является его масштабная инвариантность (фрактальность), оказывающая существенное влияние на качество связи [4-6]. Исследования трафика концентрируются вокруг статистических характеристик очередей, поскольку буферизация сообщений рассматривается как основная обеспечивающая ресурсами стратегия. Эта область исследований характеризуется в [6] следующим образом: «В 1993 г. сенсацией в области моделирования характеристик сетей стал доклад, представленный специалистами из компании ВеПСоге и Бостонского университета... Этот доклад под названием «О фрактальной природе трафика в ЕШетеЬ» по мнению некоторых специалистов, явился наиболее значительной работой по вычислительным сетям за последние десять лет. .Результаты нового взгляда на природу сетевого трафика. означают, например, что целая область проектирования компьютерных

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 3 (113) 2012 ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 3 (113) 2012

Таблица 1

Характеристики м.о. М(х) сл.в. х, реализующей сл.в. xePa(1, а)

N=106 N=1015

а M=M(x) М[х) ЛМ Ym М(х} ЛМ Ym

1,01 101 13,42 87,59 7,53 29,66 71,34 3,41

1,05 21 10,42 10,58 2,02 17,06 3,94 1,23

1,1 11 8,03 2,97 1,37 10,55 0,45 1,04

1,25 5 4,72 0,28 1,06 4,996 0,004 1,001

1,5 3 2,98 0,02 1,01 3 0,00002 1,000

Таблица 2

Характеристики дисперсии D$ сл.в. х, реализующей сл.в. xePa(1, а)

N=106 N=1015

а D(x) Щх) Ad Yd DM Ad Yd

2,01 197,0 9,94 187,1 19,8 28,26 16,99 6,97

2,05 37,2 8,33 28,9 4,5 19,78 2,09 1,88

2,1 17,4 6,78 10,6 2,6 13,41 0,58 1,29

2,25 5,8 3,95 1,8 1,5 5,58 0,05 1,03

2,5 2,2 1,94 0,3 1,1 2,22 0,002 1,002

устройств — построение буферов и управление ими — нуждается в радикальном пересмотре.. Однако среди специалистов пока нет единого мнения о том, какие математические инструменты применимы и эффективны для его исследования и прогнозирования. Их разработка должна стать следующим важным шагом в этой области».

В нашей работе [7] предложен ускоренный метод расчета буферов для очередей в сетях с фрактальным трафиком, основанный на имитационном моделировании. Эксперименты с предложенным методом показывают его высокую эффективность. Однако вопросы точного воспроизведения фрактальных распределений, строгое решение которых необходимо для дальнейшей работы над методом, оставались до сих пор в основном лишь предметом экспериментальных исследований. Настоящая статья позволяет существенно продвинуться к точному теоретическому решению этих вопросов.

Заключение. В статье получены результаты, позволяющие находить точные ответы на вопросы о погрешностях реализации распределений с тяжелыми хвостами при имитационном моделировании

сетей с очередями. Те результаты, которые получены для распределения Парето, могут использоваться в имитационном моделировании непосредственно. Что касается других распределений с тяжелыми хвостами (и не только таких), то выполненное исследование может использоваться как непосредственное руководство, методика анализа точности, с которой они реализуются в имитационном моделировании.

Библиографический список

1. Мандельброт, Б. Фрактальная геометрия природы / Б. Мандельброт. — М. : Институт компьютерных исследований, 2002. — 656 с.

2. Вишневский, В. М. Теоретические основы проектирования компьютерных сетей / В. М. Вишневский. — М. : Техносфера, 2003. — 512 с.

3. Задорожный, В. Н. Предпосылки создания фрактальной теории массового обслуживания / В. Н. Задорожный // Омский научный вестник. — 2010. — № 2(90). — С. 182— 187.

4. Шелухин, О. И. Фрактальные процессы в телекоммуникациях / О. И. Шелухин, А. М. Тенякшев, А. В. Осин ; под ред. О. И. Шелухина. — М. : Радиотехника, 2003. — 480 с.

5. Столингс, В. Современные компьютерные сети / В. Столлингс. — СПб. : Питер, 2003. — 784 с.

6. Stallings, William. Интернет и телекоммуникации / W. Stallings. [Электронный ресурс]. — URL: http://my.online. ru/it/press/cwm/1997/world.htm (Дата обращения: 13.03.2010).

7. Задорожный, В. Н. Моделирование и расчет буферов фрактальных СМО / В. Н. Задорожный, О. И. Кутузов // Имитационное моделирование. Теория и практика (ИММОД-2011) : материалы 5 й Всерос. конф. — СПб. : ЦТ СС, 2011. — Т. 1. — С. 156-161.

ЗАДОРОЖНЫЙ Владимир Николаевич, доктор технических наук, доцент (Россия), доцент кафедры «Автоматизированные системы обработки информации и управления» Омского государственного технического университета.

Адрес для переписки: e-mail: [email protected] КУТУЗОВ Олег Иванович, доктор технических наук, профессор (Россия), профессор кафедры «Автоматизированные системы обработки информации и управления» Санкт-Петербургского государственного электротехнического университета «ЛЭТИ» им. В. И. Ульянова (Ленина).

Адрес для переписки: e-mail: [email protected]

Статья поступила в редакцию 01.06.2012 г.

© В. Н. Задорожный, О. И. Кутозов

Книжная полка

51/П33

Пирумов, У. Г. Численные методы: теория и практика : учеб. пособие для бакалавров вузов по направлению подгот. «Математика. Прикладная математика» / У. Г. Пирумов [и др.] ; Моск. авиац. ин-т , Нац. исслед. ун-т. - 5-е изд., перераб. и доп. - М. : Юрайт, 2012. - 1 о=эл. опт. диск (CD-ROM). - ISBN 978-5-99161867-0.

Учебное пособие содержит все традиционные разделы, предусмотренные программой по данной дисциплине. Материал дается по единой схеме, включающей в себя постановку задачи, описание алгоритма решения, детально разобранные типовые примеры и тщательно подобранный комплекс задач для самостоятельного решения.