2007-uch-posob-vargapetyan-1
.pdf
Математическое ожидание БСВ (первый начальный момент)
μ = E{ α} =1/2;
дисперсия (второй центральный момент)
σ2 = D{α} = E{(α −μ)2} =1/12.
Наряду с простейшим экспериментом будем рассматривать состав ной случайный эксперимент, получающийся в результате r кратно го (r ≥1) повторения независимо друг от друга простейших экспери ментов. Результатом составного случайного эксперимента является последовательность из r независимых БСВ α1, …, αr таких, что
α |
i |
= α |
(ω) = ω , i = |
1, r; |
ω=(ω ,..., ω ) Ωr , |
|
|
i |
i |
1 |
r |
||
где ωi — координата точки, брошенной наудачу в [0, 1) в i м простей шем эксперименте.
Совместная плотность распределения вероятностей α1, …, αr
|
|
|
|
|
|
|
|
|
|
|
|
|
|
[0,1), i =1, r; |
|||
p |
|
|
|
1, x |
||||
|
(x ,..., x ) = |
i |
|
|
|
|||
α1 |
,...,αr |
1 |
r |
0 в противном случае. |
||||
|
|
|
|
|
|
|
|
|
Согласно второму принципу моделирования случайных элемен тов, любой СЭ Θ представляется для некоторого натурального r в виде функции f(.) от r независимых БСВ:
Θ = f(α1, …, αr).
Таким образом, задача моделирования произвольного СЭ Θ* раз бивается на две подзадачи:
1)моделирование на ЭВМ независимых БСВ α1, …, αr;
2)нахождение функции f(.) такой, чтобы СЭ Θобладал требуемы ми вероятностным законом распределения и числовыми характерис тиками.
Поэтому моделирующий алгоритм состоит из двух блоков
(рис. 3.8):
Для имитации одного и того же СЭ Θ* может быть предложено несколько вариантов функциональных преобразований. Обычно пред
Б1 |
α1, …, αr |
Б2 |
Θ |
|
|
||
|
|
|
|
Рис. 3.8. Моделирующий алгоритм БСВ: Б1 — блок моделирования БСВ (об( щий для всех Θ); Б2 — блок функционального преобразования f(.) БСВ (различный для различных законов распределения вероятностей)
61
почтение отдается варианту f(.), требующему меньших вычислитель ных затрат; для этого применяется понятие коэффициента исполь зования БСВ.
Коэффициентом использования БСВ K назовем величину, обрат ную числу r базовых случайных величин, используемых для модели рования одной реализации случайного элемента Θ*:
K =1
r, 0 < K ≤1.
Величина K является мерой вычислительных затрат на модели рование Θ*. Чем меньше K, тем больше затраты. Целесообразно вы бирать такую функцию f(.), для которой K принимает наибольшее значение.
Очевидно, чтобы моделировать на ЭВМ случайные элементы с за данным вероятностным законом распределения, необходимо уметь моделировать БСВ. БСВ αявляется абсолютно непрерывной случай ной величиной. Однако на ЭВМ приходится иметь дело с дискретны( ми случайными величинами (ДСВ). Поэтому моделирование БСВ основано на аппроксимации непрерывной случайной величины αДСВ α′. Опишем способ построения ДСВ α′.
Рассмотрим случай, когда представление целых неотрицательных чисел на ЭВМ осуществляется с помощью k двоичных разрядов (би тов). Тогда С = {0, 1, 2k – 1} — множество 2k неотрицательных целых чисел, представимых в ЭВМ. Определим на (Ω, Φ, Р) ДСВ β = β(ω) следующим образом:
0, 0 ≤ ω≤2−k; |
|
|
|
|||||
|
|
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
≤ ω< (i +1)2−k, i C; |
(3.18) |
||||
β =β(ω) = i =i2−k |
||||||||
|
|
|
|
|
|
|
|
|
1 |
|
|
|
−k |
|
−k |
|
|
|
k |
−1,1 |
−2 |
k |
=1. |
|
||
2 |
|
≤ ω<2 2 |
|
|
||||
|
|
|
|
|
|
|
|
|
Построение βпроиллюстрируем с помощью рис. 3.9. Разобьем про межуток [0, 1) на 2k отрезков одинаковой длины 2–k; для ω, попадаю щих в промежуток [i2–k, (i + 1)2–k), полагаем β =β(ω) =i, i C.
ω
0 2–k 2•2–k … i2–k (i + 1)2–k … 1–2–k 1
Рис. 3.9. Построение β
62
По построению распределение ДСВ β является равномерным на множестве С, т. е. все значения βравновероятны, действительно:
P{β(ω) =i} = P{ω:ω i2−k, (i +1)2−k )} =2−k.
Теперь перейдем от СВ βк искомой ДСВ α′:
α′= α′(ω) =β(ω)/2k =β(ω)2−k.
Согласно (3.19):
0, 0 ≤ ω≤2−k, (β=0);
1
α′= −k −k ≤ ω<( + ) −k β =
i2 , i2 i 1 2 , ( i), i C;
1
1−2k, 1−2−k ≤ ω<1, (β =2k −1),
(3.19)
(3.20)
т. е. от целочисленной ДСВ β мы перешли к ДСВ α′ со значениями в [0, 1).
Очевидно, все возможные значения α′определяются множеством С′ = {0, 2–k, …, 1–2–k} и являются равновероятными:
P{α′(ω) =i2−k} = P{β(ω) =i} =2−k, i C,
т. е. закон распределения α′является равномерным на C′.
Точность аппроксимации αс помощью α′устанавливается с помо щью леммы.
Лемма 3.1. Для СВ α= α(ω) и α′ = α′(ω), определенных на (Ω, Φ, Р) и имеющих вид (3.16), (3.17) соответственно, равномерное уклоне ние удовлетворяет выражению
sup | α′(ω) −α(ω) |<2−k. |
(3.21) |
ω Ω |
|
Доказательство. Разобьем Ω= [0, 1] на 2k промежутков согласно рис. 3.9. Пусть ω i2−k, (i +1)2−k ), i C. Тогда, согласно (3.16), спра ведливо представление
α=α(ω) =ω=i2−k +δ, 0 ≤δ<2−k, i C,
асогласно (3.20):
α′=α′(ω) =i2−k.
Поэтому α−α′ = δ<2−k, ω Ω. Отсюда заключаем справедливость
(3.21).
63
Fα(x), Fα(x)
1 
1 – ε
ε
|
|
|
0 |
ε ... 1 – ε |
1 |
x |
|
|||
Рис. 3.10. Функция распределения Fα(x) |
|
|
|
|
||||||
Таблица 3.4. Соотношения между дисперсиями величин α и α′ |
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
k |
|
2 |
|
3 |
|
5 |
|
10 |
15 |
|
|
|
|
|
|
|
|
|
|
|
|
D{α } |
1.290 |
|
1.140 |
|
1.030 |
|
1.001 |
1.00 |
|
|
′ |
1/2 |
|
|
|
|||||
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
D{α} |
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Из (3.21) следует, что если k → ∞, то последовательность α′k → α равномерна по ω Ω. Таким образом, случайная величина α′ = α′(ω) является аппроксимацией для БСВ1 α=α(ω); α′называется в связи с этим квазиравномерной случайной величиной. Ее функция распре деления Fα′ (x) показана на рис. 3.10 и аппроксимирует Fα (x) с точ ностью ε = 2–k.
Между математическими ожиданиями величин α, α′справедливо соотношение
E{α′} =0,5(1−2−k ) → E{α} =0,5,
представленное в табл. 3.4.
При достаточно больших значениях k (например, для ПЭВМ IBM PC AT 486, Pentium k = 31, ε= 0.5•10–9) величины αи α′отождеств ляют.
3.4.3. Методы построения программных генераторов
Существует ряд методов и алгоритмов, используемых при постро ении генераторов БСВ [4, 11]. Наиболее широкое распространение получили генераторы, построенные на использовании различных
1 Напомним, что аппроксимация непрерывных случайных величин диск ретными широко используется в теории вероятностей при построении интегра ла Лебега.
64
модификаций метода вычетов (мультипликативный конгруэнтный (совпадающий) метод), который в разных модификациях использует ся в разных редакциях GPSS/H. В конце параграфа будут приведены параметры генераторов GPSS/H. Рассмотрим этот метод подробнее.
Любое следующее число образуется на основании выражения
xi+1 = axi + c mod m i = 1, 2, …, |
(3.22) |
||
где a, m — положительные целые числа; c — целое число. |
|
||
Исторически этот генератор строился из принципа |
|
||
x |
n |
= xn mod m. |
(3.23) |
|
|
|
|
Чтобы не возникало вычислительных проблем, (3.23) можно заме нить
xn = xxn–1 mod m. |
(3.24) |
При с = 0 как раз и получается мультипликативный конгруэнт ный генератор и генерируемые числа Ui = xi / m. Методы проверки случайности генерируемых чисел будут рассмотрены в п. 3.5.3. Вы ражение (3.24) основывается на значении m = 2k, где k — целое чис ло, обычно равное длине машинного слова. Любые генерируемые БСВ характеризуются периодом повторения или циклом генератора d = 2k–2 при нечетном x0 (примем это без доказательства).
В качестве примера рассмотрим генератор xn = 3xn–1 mod 32 (25). При с = 0 [см. (3.22)] и а = ±3 (mod 8) период этого генератора d = 85 – 2 = 8. Приняв x0 = 1, найдем числовую последовательность,
создаваемую генератором:
x0 = 1 |
00001 |
x1 = 3 |
00011 |
x2 = 9 |
01001 |
x3 = 27 |
11011 |
x4 =81 mod 32 (81– 2*32) =17 |
10001 |
x5 =51 mod 32 = 19 |
10011 |
x6 =57 mod 32 = 25 |
11001 |
x 7 = 75 mod 32 =11 |
11001 |
x8 = 33 mod 32 = 1 |
00001 |
Последовательность повторилась, т. е. период равен 8, как и было подсчитано, в случае четного x0 период сокращается вдвое!
Современные ЭВМ (Pentium 1 и выше) сняли все ограничения по быстродействию и объемам памяти, необходимым для генерации ка чественных БСВ. Потому при выборе генератора БСВ надо руковод ствоваться следующими правилами:
1) период повторения должен быть более миллиона;
65
2)БСВ, производимые генератором, должны отвечать всем тестам на случайность;
3)стараться пользоваться стандартными генераторами типа RANDU (URN01 — URN42), встроенными непосредственно в ЯИМ.
Рассмотрим более подробно схему повторения БСВ. На интервале
(0, m) должны быть такие числа r и s (0 ≤ r < s < m) , для которых xr =xs, а следовательно, иxr+l = xs+l для всех положительных l. Если s наименьшее положительное число, для которого xr = xs, тогда выра жение (3.22) состоит из начального сегмента x0, x1, …, xr–1, предше ствующего второму сегменту xr, xr+1, …, xr = xs,который затем посто янно повторяется. Тогда целое число d = s – r называется периодом повторения, последовательностью или циклом генератора. После довательность имеет максимальный период, если d = m. Проиллюст рируем зависимость периода d от a, c, m, x0 на двух простых приме рах.
Пример 1. Пусть x0 = 1, a = 5, c = 0, m = 10, тогда определяем последо вательность x1= 9 mod 10 = 9, x2 = 81 mod 10 = 1, x3 = 9 mod 10 =9, …. Последовательность имеет вид 1, 9, 1, 9, … . В данном случае r = 0, s = 2 и период d = s – r = 2.
Пример 2. Пусть x0 = 1, a = 5, c = 0, m = 10, тогда определяем последо вательность x1= 5 mod 10 = 5, x2 = 25 mod 10 = 5, …. Последовательность имеет вид 1, 5, 5, …. В данном случае r = 1, s = 2, d = 1. У данного генератора начальный сегмент 1 предваряет последовательное число 5.
В табл. 3.5 приведены данные генератора для различных значе ний x0, a для уравнения (3.23) — при с = 0 и m = 10.
Таблица 3.5. Данные генератора
x0 |
|
|
|
|
a |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
||
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
1 |
1 |
1,2,4,8,6, |
1,3,9,7, |
1,4,6, |
1,5, |
1,6, |
1,7,9,3, |
1,8,4,2,6, |
1,9, |
|
|
|
|
|
|
|
|
|
|
|
|
2 |
2 |
2,4,8,6, |
2,6,8,4, |
2,8, |
2,0, |
2, |
2,4,8,6, |
2,6,8,4, |
2,8, |
|
|
|
|
|
|
|
|
|
|
|
|
3 |
3 |
3,6,2,4,8, |
3,9,7,1, |
3,2,8, |
3,5, |
3,8, |
3,1,7,9, |
3,4,2,6,8, |
3,7, |
|
|
|
|
|
|
|
|
|
|
|
|
4 |
4 |
4,8,6,2, |
4,2,6,8, |
4,6, |
4,0, |
1,4, |
4,8,6,2, |
4,2,6,8, |
4,6, |
|
|
|
|
|
|
|
|
|
|
|
|
5 |
5 |
5,0, |
5, |
5,0, |
5, |
5,0, |
5, |
5,0, |
5, |
|
|
|
|
|
|
|
|
|
|
|
|
6 |
6 |
6,2,4,8, |
6,8,4,2, |
6,4, |
6,0, |
6, |
6,2,4,8, |
6,8,4,2, |
6,4, |
|
|
|
|
|
|
|
|
|
|
|
|
7 |
7 |
7,4,8,6,2, |
7,1,3,9, |
7,8,2, |
7,5, |
7,2, |
7,9,3,1, |
7,6,8,4,2, |
7,3, |
|
|
|
|
|
|
|
|
|
|
|
|
8 |
8 |
8,6,2,4, |
8,4,2,6, |
8,2, |
8,0, |
1,8, |
8,6,2,4, |
8,4,2,6, |
8,2, |
|
|
|
|
|
|
|
|
|
|
|
|
9 |
9 |
9,8,6,2,4, |
9,7,1,3, |
9,6,4, |
9,5, |
9,4, |
9,3,1,7, |
9,2,6,8,4, |
9,1, |
|
|
|
|
|
|
|
|
|
|
|
66
По данным табл. 3.5, наибольшее значение периода d равно 4, а плохой выбор a и x0 значительно уменьшает этот период, в таблице подчеркнуты повторяющиеся сегменты.
Для выбора генератора, обладающего наибольшим периодом, рас смотрим несколько лемм.
Лемма 3.2. Если a и m — простые числа, всегда есть положитель ное целое d, для которого x0 = xd. Для простых чисел наибольшим общим делителем (НОД) является единица.
Лемма 3.3. Если НОД (a, m) = 1, то тогда период последователь ности (3.22) является наименьшим целым положительным числом d, для которого
σd ((a – 1) x0 + c) = 0 (mod m), |
(3.25) |
где σd = 1 + a + a2 + …+ a d–1.
Необходимо отметить, что эти леммы действительны только для случая, когда НОД равен 1, поэтому для гарантии надо выбирать только такой генератор, у которого НОД равен 1.
Стремление увеличивать период является необходимым, но не достаточным условием. При увеличении периода может произойти ухудшение качества БСВ, поэтому их необходимо проверять на слу чайность, что и рассмотрено в § 3.5.
В заключение параграфа рассмотрим тип генератора БСВ, приме няемый в GPSS/H. В первой версии языка использован генератор URN 27, работающий на сдвиговом принципе с обратной связью и использующий выражение x31 + x3 + 1. Это выражение комбинирует ся с последовательностью xi = 69069 xi–1 mod 232. Сдвиг вправо на 15 разрядов и влево на 17 обеспечивает получение 32 разрядного числа yi. Эти xi и yi складываются по правилам алгебры логики и создают последовательность zi, так, например:
xi = 01001…01010 yi = 00101…10110 zi = 01100…11100
Этот генератор не проходил по некоторым тестам случайности, по этому в более поздних версиях использован усовершенствованный ге нератор, имеющий в качестве минимального 32 разрядного числа (по ложительный ноль и 31 значащий разряд) минимальное значение 1 и максимальное 2 млрд 147 млн 483 тыс. 646 (231 – 2), а выражение принимает вид xi+1 = 742938285 xi mod 2147483647, при этом x0 = 266301881. Таким образом, период повторения этого генератора весьма велик (более 2 млрд), что вполне достаточно для моделирова ния сложных систем. Независимость потоков БСВ достигается тем,
67
что потоки заявок и обслуживания берутся с разных генераторов. В ранних версиях использовались 8 встроенных генераторов, старту ющих с одинаковых начальных позиций, что делало потоки идентич ными, а задачу сбора статистики достаточно сложной. В последних версиях статистическая независимость достигается практически не ограниченным числом потоков и отличием каждой последующей вы борки на 100 000. Такое отличие потоков практически исключает воз можность их корреляции. Более тонкие моменты исследования регу лярности этих процессов в пособии не рассматриваются. Кроме того, регулировка сдвига начальных значений легко осуществляется с по мощью оператора управления (ОУ) RMULT и ОБ BRMULT, которые позволяют сдвигать начальные значения на один миллион и более! Следует подчеркнуть некоторые методические особенности использо вания генераторов БСВ. Разброс времени A±B в ОБ GENERATE и ADVANCE всегда берется с генератора RN1, с этого же генератора бе рется случайное значение ОБ TRANSFER в статистическом виде.
§ 3.5. ОЦЕНКА КАЧЕСТВА БАЗОВЫХ СЛУЧАЙНЫХ ВЕЛИЧИН, ПОЛУЧАЕМЫХ ОТ ПРОГРАММНЫХ ГЕНЕРАТОРОВ
3.5.1. Общие представления
Рассмотрим кратко вопросы оценки качества БСВ, т. е. соответ ствие последовательности БСВ равномерности распределения в ин тервале (0,1). Такую оценку позволяют осуществлять статистиче ские тесты. Поскольку причин отклонения БСВ от случайности мно го, то и тестов проверки также много [15]. Идеальным генератором считается тот, который проходит все статистические тесты. Каждый из таких тестов ориентирован на определение частной причины от клонения последовательности от независимости и случайности. Бу дем считать, что существует k возможных причин отклонений, и при выполнении гипотезы случайности вероятность появления i(й при чины равняется πi, i = 1, 2, …, k. Предположим, что рассматривается достаточно большая последовательность, в которой по разным при чинам возникло n несоответствий по всем k причинам. Количество несоответствий по каждой группе обозначим как fi, тогда f1 + …+ fk = n; так как fi имеет вероятность πi, то число несоответствий равно nπ1 = ei . Измерение различия между fi и ei определяется из выражения
|
|
|
|
|
|
|
k |
|
|
D = (f – e |
)2/ e |
1 |
+ … + (f |
– e |
)2 / e |
k |
= ∑(f −e )2 |
/e . |
(3.26) |
1 1 |
|
k |
k |
|
i i |
i |
|
||
|
|
|
|
|
|
|
i=1 |
|
|
68
Таблица 3.6. Значения с для различного числа степеней свободы γ
|
1 |
2 |
3 |
4 |
5 |
10 |
20 |
50 |
100 |
|
|
|
|
|
|
|
|
|
|
с |
3.841 |
5.992 |
7.814 |
9.487 |
11.071 |
18.307 |
31.411 |
67.504 |
124.342 |
|
|
|
|
|
|
|
|
|
|
При истинной случайности разность между fi и ei должна стре миться к нулю, а D должно быть достаточно малым, следовательно, гипотеза о случайности отвергается при большом значении D. При уровне значимости 0.05 вероятность отвергнуть истинную гипотезу о случайности равна
PH0 (D > c) = 0,05. |
(3.27) |
Оценить точное распределение D достаточно сложно, однако при n > 5 D приблизительно распределено по хи квадрат распределению с числом степеней свободы γ= k – 1. Для большей уверенности общее число несоответствий обычно принимают равным 10. В табл. 3.6 при ведены значения c для различных значений числа степеней свободы с тремя значащими цифрами после точки, при которых справедливо равенство (3.27).
Тесты позволяют исключить процедуру статистического оценива ния, при которой из за ошибок 1 го и 2 го рода оценка вероятности принятия правильной гипотезы равномерности распределения ока зывается достаточно сложной. Ниже кратко рассмотрим два типа тестов: тесты на равномерность БСВ и тесты проверки качества са мих генераторов.
3.5.2.Тесты оценки качества БСВ
1.Тест равномерности распределения (хи!квадрат!тест).
Второго названия этого теста следует избегать, чтобы не возникло
путаницы из за совпадения названий при рассмотрении распределе ния D. При этом тесте интервал (0, 1) делится на сто равных и не пере секающихся интервалов длиной 0.01, вероятности появления несо ответствий в каждой подвыборке из n равны e1 = e2 = … = e100 = n/100.
Затем проверяется равенство вероятностей в каждом интер вале.
2. Тест замещения (сбора купонов).
Имеем последовательность случайных чисел U1, U2, …, произве дем замену чисел на целочисленные значения 1, 2, …, M по правилу:
Заменим Ui на 2, если 1/M ≤ Ui < 2/M
69
Заменим Ui на 1, если 0.0 ≤ Ui < 1 /M
.
.
.
Заменим Ui на M, если (M – 1)/M ≤ Ui ≤ 1.0
Эти новые целые числа от 1 до М должны быть случайными, если Ui действительно случайные числа. Рассмотрим последовательно всю полученную совокупность и обозначим через Q новые числа, полу ченные после М, например М + 1, М + 2, …,. При соблюдении слу чайности вероятности Р (Q = M), P (Q = M + 1), P (Q = M + 2), … известны. Если вычислить такие значения Q повторно, то тогда тест 2 может быть применен для проверки случайности начальных
чисел Ui.
3. Тест расхождения.
Положим, что αи βлежат между 0.0 и 1.0 при условии, что α< β. Имеющуюся последовательность чисел U1, U2, …, будем рассматри вать с позиции их нахождения: либо в интервале [α, β] включитель но границы, либо вне его. Будем считать, что каждое число, лежащее в интервале, равно 1, а вне его равно 0. Эта операция преобразует последовательность в набор 0 и 1. Положим, что p = β– α, тогда, если последовательность чисел случайна, то вероятность того, что j(й 0 появится после 1, прежде чем появится следующая 1, равняется
p |
j |
= p (1 – p) j, j = 0, 1, 2,… |
(3.28) |
|
|
|
После определения, сколько раз появится каждая категория, мож но применить тест 1, который оценивает случайность в соответствии с (3.28). Этот тест введен Кендаллом и Бабингтон Смитом и предус матривает следующие правила для выбора αи β:
вариант 1 α = 0.0, β = 0.5 (ниже среднего); вариант 2 α = 0.5, β = 1.0 (выше среднего);
вариант 3 α = 0.333, β = 0.667 (нахождение в середине).
4. Тест перестановки.
Считается, что полученная последовательность случайных чисел U1, U2 , …, содержит наборы T чисел. Каждый такой набор содержит T! возможных перестановок, которые можно рассматривать как T чисел, расположенных в порядке убывания от наибольшего к наи меньшему. При случайности вероятность появления каждого набо ра T! равна 1/T!. При T! категориях можно определить, сколько кор тежей образовано из соответствующих наборов T и сколько из них попало в соответствующую T!(категорию. Затем применяется тест 1, для этого теста рекомендованы значения T = 3, 4, 5.
70
