Имитационное моделирование экономических процессов

пример 1.1. Предположим, что нужно проверить данные, полу­ ченные (или наблюдаемые) при использовании метода Монте-Карло и приведенные в табл. 1.1, на их соответствие распределению Пуас­ сона.

Таблица 1.1

Данные для проверки гипотезы по критерию Колмогорова-Смирнова

Число Наблю­ собы­ даемая тий частота

1

0 315

1 142

240

39

42

51

Наблю-,	Теоре­	Наблю­	Теорети­
даемая	тическая	даемое	ческое
вероят­	вероят­	распреде­	распреде­
ность	ность	ление	ление
2	3	4	5
0,619	0,571	0,619	0,571
0,279	0,319	0,898	0,890
0,078	0,089	0,976	0,979
0,018	0,017	0,994	0,996
0,004	0,003	0,998	0,999
0,002	0,001	1,000	1,000

Абсо­

лютная

разность

6

0,048

0,008

0,003

0,002

0,001

0,000

Эти данные имеют следующий смысл. На отрезке времени наблюдаем случайные собьггия, число Kotopbix равно х. Если это распределение Пуассона, то вероятность Р{х = п} того, что х = п, где п - заданное число, равна

Р{х = и} =Ге-^

где е= 2,71828; Х- положительная константа, которая одновременно является и математи­

ческим ожиданием, и дисперсией.

Предположим, что Я, = 0,5577. Сформулируем гипотезу в сле­ дующем виде: Не имеется существенных различий между наблю­ даемыми данными во время эксперимента и теми данными, которые должны получаться из распределения Пуассона расчетным путем со средним значением 0,5577 и JV= 509.

Рассмотрим подробнее табл. 1.1. В ней строки с номерами О,..., 5 соответствуют числу п наблюдаемых событий, и=0, ... , 5, столбец 1 содержит наблюдаемую частоту появления и событий, а столбец 2 - наблюдаемую вероятность появления и событий. В столбце 3 пред­ ставлены рассчитанные по формуле значения вероятности появления п событий.

21

Прежде всего нужно получить два интегральных распределения (фактически приближения функций распределения). Сначала сдела­ ем это для наблюдаемых данных: с помощью столбца 2 получим столбец 4. Затем - для теоретических данных: с помощью столбца 3 получим столбец 5.

После этих вычислений найдем абсолютные разности для всех групп значений случайной величины и с помощью столбцов 4 и 5 получим столбец 6. В последнем столбце наибольшая абсолютная разность 0,048 получается в группе, соответствующей нулевому числу событий.

Далее необходимо найти так называемое критическое значение Door для проверки принятой гипотезы. Таблица критических чисел многократно переиздавалась. Критические числа в виде, удобном для вьшолняемой проверки, приведены в табл. 1.2. Абсолютную разность 0,048 необходимо сравнить с критическим значением, най­ денным по табл. 1.2.

					Т а б л и ц а 1.2
	Критические числа Колмогорова-Смирнова
Степень		Проверка		Проверка
свободы	единичной выборки *			двух выборок **
N	А,10	^0,05	^0,01	Dofls	^0,01
1	0,950	0,975	0,995	-	-
2	0,776	0,842	0,929	-	-
3	0,642	0,708	0,828	-	-
4	0,564	0,624	0,733	1,000	1,000
5	0,510	0,565	0,669	1,000	1,000
6	0,470	0,521	0,618	0,833	1,000
7	0,438	0,486	0,577	0,857	0,857
8	0,411	0.457	0,543	0,750	0,875
9	0,388	0,432	0,514	0,668	0,778
10	0,368	0.410	0,490	0.700	0,800
11	0,352	0,391	0,468	0,636	0,727
12	0,338	0.375	0,450	0,583	0,667
13	0,325	0.361	0,433	0,338	0,692
14	0,314	0,349	0,418	0,571	0,643

22

						Продолжение
"Степень		Проверка		Проверка
свободы	единичной выборки *			двух выборок
N	^0,10	^0,05	•Do.oi	^^0,05		•Do.Ol
15	0,304	0,338	0,404	0,533		0,600
16	0,295	0,328	0,392	0,500		0,625
17	0,286	0,318	0,381	0,471		0,588
18	0,278	0,309	0,371	0,500		0,556
19	0,272	0,301	0,363	0,474		0,526
20	0,264	0,294	0,356	0,450		0,550
25	0,240	0,270	0,320	0,400		0,480
30	0,220	0,240	0,290	0,370		0,430
35	0,210	0,230	0,270	0,340		0,390
Более 35	1,22	1,36	1,63	1,3б1^1^^2	1.63	1Ш^^
		7W
				V NiNz

* Применяется для оценки степени близости выборочных значений к тео­ ретическому распределению.

N - объем выборки.

** Применяется для определения принадлежности двух выборок объемами

Ni и JVj одному и тому же распределению. При малых размерах выборки

N^Ny'Ni.

При и=509 и значении индекса критического числа Da а = 0,05 получается критическое значение

Z ) ^ = b 3 6 ^ U 6 ^ ^ J [ ^ ^'"^ 7N 7509 22.56 '

Поскольку наибольшая разность 0,048 < Z5extr. то не отказьюаемся от гипотезы о том, что экспериментальное распределение - пуассоновское. Проверка статистических гипотез о соответствии «событий явлению» и «явлений поведению» дает математический инструмент для оценки «рискованного поведения» исследуемого процесса.

23

1.3

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

ЭКОНОМИЧЕСКИХ ПРОЦЕССОВ

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

Равномерное распределение на интервале (0,1). В литературе приводились описания разных датчиков случайных величин для получения последовательностей чисел, распределенных по какомуто случайному закону. Основная проблема заключалась в программ­ ной реализации равномерного распределений на интервале (О, 1). Существуют различные методы получения такого равномерного распределения. Остановимся на программном генераторе, наиболее подходящем для компьютеров с 32-разрядным словом. Период по­ следовательности, получаемой с помощью такого генератора, на несколько порядков превосходит период, получаемый по методу, рассмотренному в разд. 1.2.

Рассмотрим текст функции на языке C++, реализующей этот ге­ нератор и возвращающей величину, равномерно распределенную на интервале (0,1). Данная функция реализует одну из разновидностей мультипликативного конгруентного метода. Генератор предназначен для применения в системе имитационного моделирования, позво­ ляющей параллельное моделирование сложной сети взаимодейст­ вующих процессов, причем каждый процесс может иметь свой дат­ чик псевдослучайных величин. Поэтому в качестве глобальной пе­ ременной рассматривается указатель к - адрес управляющий струк­ туры такого процесса, имеющего номер next. Ниже следует текст программы.

Все математические приемы и программные датчики псевдослучайных величин, рассмотренные в данном разделе, имеются в составе системы Pilgrim.

24

•float rundum(void)

{

s t r u c t kcb *k;

к= addr[next];

iy = k->bx * 1220703125;

if	(iy	< 0)	+ 1073741824;
rn	iy = iy + 1073741824
	= iy	* 0.4656613e-9;
k-»bx =		iy;

r e t u r n ( r n ) ;

}

Внутри этой структуры расположено 32-разрядное поле Ьх, ко­ торое адресуется как к->Ьх. Перед моделированием управляющая программа должна занести в это поле длинное нечетное число, на­ пример k->bx=(long)2098765431. Следует отметить, что для каждого процесса нужно задать свое число. Этого можно достигнуть, напри­ мер, применяя для каждого следующего инициируемого процесса начальное значение предьщущего, уменьшенное на 2. Если мы хо­ тим резко повысить эффективность работы этой функции и еще сильнее приблизить псевдослучайную последовательность к случай­ ной, то можно в качестве начального значения использовать бито­ вую последовательность таймера ЭВМ: k->bx=(long)time(NULL). Далее это число проверяется на нечетность, и если оно четное, то его необходимо увеличить на 1.

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

•обратных функций;

•комбинирования величин, распределенных по другим законам, например по равномерному на интервале (0,1).

Равномерное распределение на произвольном интервале.

Рассмотрим важное и очень простое равномерное распределение на интервале (m-s, m+s). Плотность вероятностей этого распределения описьшается следующей формулой:

	О	,	t<m-s	;
Pit)^	1	,	m-s<t<m	+ s
	2s		t>m + s
	О	,

где m - математическое ожидание;

i - максимальное отклонение от математического ожидания.

25

Такое распределение используется, если об интервалах времени известно только то, что они имеют максимальный разброс, и ничего не известно о распределениях вероятностей этих интервалов.

Программная функция на языке C++ приведена ниже (два ее входных параметра не нуждаются в комментариях):

float unifrm( float m, -float s) f l o at a,b;

a= m-s;

b= m+s;

return(a+ (b-a) *runduin() ) ;

}

Равномерное распределение можно использовать при расчетах по сетевым графикам работ, в том числе при работе по методу PERT. Это распределение можно применять и при расчетах основньк дли­ тельностей и времен в военном деле (времени выдвижения воинской части или ее подразделения на исходный рубеж, времени марша, времени подготовки рубежа обороны и др.).

Формула для определения дисперсии получается после получе­ ния второго момента:

nrj^-, Aim+s)-{m-s)f		_s^
^ ^	12	3 •

Дискретное распределение (общий случай). Предположим, что известны частоты а/ выбора из N объектов на определенном интер­ вале времени, »-1,..., N. Пример таких частот для N=1 представлен в табл. 1.3. Первая строка таблицы - это номер объекта, а вторая - частота его выбора. Требуется разработать программную функцию, которая должна возвращать значение номера объекта в соответствии с этими частотами.

						Т а б л и ц а 1.3
		Значения обратных функций для получения
		дискретного распределения
1	1	2	-3	4	5	6	7
«/	130	5	И	44	32	2	67
Рг	0,447	0,017	0,038	0,151	0,110	0,007	0,230
Y*	0,447	0,464	0.502	0,653	0,763	0,770	1,000

26

Воспользуемся методом обратных функций. Сначала найдем сумму всех частот:

N

1=1

В нашем случае получаем Л=291. После этого построим таблицу нормированных значений Р,=а,/Л (третья строка табл. 1.3). Далее рассчитаем значения дискретной функции у, по формуле

к=1

Значения у, находятся в четвертой строке табл. 1.3.

Построим график дискретной функции у, (рис. 1.1). Далее вос­ пользуемся рассмотренной выше программой получения случайных величин, распределенных равномерно на отрезке (0,1), и каждый раз будем получать случайную величину/>,. Условимся, что Уо = О . По­ сле этого выбор объекта с номером / осуществляется при вьшолнении соотношения

у,_1<А^У/-

У' i 1

1

Рл

						Р,
	Pt -			^-	•Р4
				Р\
		* •	t^P2	'
	о	Pi-.
		1	2	3	4	i	Л^

Номер объекта

Рис. 1.1. Пример применения метода офатных функций

27

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

Непрерывная случайная величина t имеет нормальное распреде­ ление вероятностей с параметрами ти и а > О, если ее плотность ве­ роятностей имеет вид (рис. 1.2):

1

Р(0 = -г>/2я ехр

где т - математическое ожидание M[i]; о - среднеквадратичное отклонение.

Рис. 1.2. Вид нормального распределения

Причем если D[t] - это дисперсия, то о = VДО •

Рассмотрим временные диаграммы на рис. 1.3. Предположим, что какой-то случайный процесс состоит из последовательности и

28

элементарных независимых процессов. Длительность каждого эле­ ментарного процесса Г/ - это случайная величина, распределенная по неизвестному закону с математическим ожиданием /, и дисперсией

aj. Допустим, что это непрерывное распределение (допущение,

справедливое для нормальной жизни общества, где существует при­ родная инертность), причем третий момент должен иметь ограниче­ ние по абсолютной величине (это также, естественно, реальное ус­ ловие). Справедливо соотношение

Пп} = 1г,

/=1

t„

О

п - случайная величина, являющаяся суммой п неза- Т{п}= Х^- висимых случайных величин, распределенных ло i-l неизвестному закону и имеющих конечный третий

абсолютный момент.

Если сделать предельный переход п->оо, то распределение вели­ чины t=T{n} стремится к нормальному с математическим ожида­ нием М[(1, равным сумме математических ожиданий элементарных случайных отрезков Г/, i=l,..., п, и дисперсией D[t], равной сумме дисперсий этих отрезков:

/=1 /=1

Рис. 13. Интерпрета'щя нормального распределения

Теорема 1 (без доказательства). Если сделать предельный пере­ ход и устремить И-+00 , то распределение случайной величины t=T{n} устремится к нормальному с математическим ожиданием Щ(\ и дисперсией D[i\, определяемыми из следующих соотношений:

M[t\=t~t.y^D{t\=tah

1=1 1=1

В различных математических постановках центральная предель­ ная теорема рассматривается в научной литературе по теории веро­ ятностей и математической статистике.

29

практический смысл этой теоремы очень прост. Любые сложные работы на объектах экономики (ввод информации из документа в компьютер, проведение переговоров, ремонт оборудования и др.) состоят из многих коротких последовательных элементарных со­ ставляющих работ. Причем количество этих составляющих работ иногда настолько велико, что требования в приведенной вьнпе тео­ реме о независимости и- одинаковом распределении становятся из­ лишними. Поэтому при оценках трудозатрат всегда справедливо предположение о том, что их продолжительность - это случайная величина, которая распределена по нормальному закону.

Ниже приведена функция на C++, возвращающая нормально распределенную случайную величину:

float normal(float m, float s)

{

float a; int i ; a=0.0;

for (i=0; i<12; i++) a += rundijmiO ;

return(m+{a-6.0)* s) ;

}

Работа функции основана на применении центральной предель­ ной теоремы. При получении одного числа используются 12 равно­ мерно распределенных на отрезке (О, 1) величин, которые суммиру­ ются. Последовательность чисел, распределенных равномерно на отрезке (0,1), имеет математическое ожидание 1/2 и дисперсию 1/12. С учетом центральной предельной теоремы сумма таких 12 чисел имеет математическое ожидание 6 и дисперсию 1. После суммиро­ вания выполняются необходимые действия для обеспечения пара­ метров нормального распределения: математического ожидания т и дисперсии о^.

Входными параметрами этой функции являются:

•m - математическое ожидание M[t];

•S - среднеквадратичное отклонение а = 7 д ? ] .

Преимущество этой функции - высокое быстродействие. Недос­ татком является игнорирование «хвостов» нормального распределе­ ния, которые могут уходить в обе стороны от величины т на рас­ стояние, превышающее 6а. Поэтому при проведении особо точных экспериментов применяются другие - более точные (но более мед­ ленные) функции. В современных системах имитационного моде-

30