Инструментарии статистического моделирования систем

Сущность метода статистического моделирования заключается в том, чтобы с помощью специального моделирующего алгоритма воспроизвести на ЭВМ поведение системы и процесс функционирования ее элементов с учетом воздействий случайных факторов внутреннего и/или внешнего происхождения. В результате применения этого весьма мощного аналитического метода исследования и моделирования систем получается серия частных значений искомых функций и величин, статистическая обработка которых с помощью методов теории вероятностей и математической статистики позволяет получить информацию о поведении системы в произвольные моменты времени. Если в процессе экспериментирования с моделью накоплены данные достаточно большого объема (т. е. количество выборки достаточно велико), то результаты их обработки приобретают статистическую устойчивость и с достаточной точностью могут быть приняты в качестве надежных оценок искомых характеристик системы и ее поведения с учетом взаимодействия с внешней средой.

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

а)Теорема Бернулли. Если некоторое событие A в N независимых экспериментах осуществляется с вероятностью P, то относительная частота появления этого события m/N сходится по вероятности к P, т.е. для любого > 0 имеет место

lim P_r{|m/N – P| ³ e} = 0, (4.1)

N® ¥

где m – число положительных исходов в N экспериментах.

б)Теорема Пуассона. Если в N независимых испытаниях события A появляется с вероятностью в каждом i - м испытании, то относительная частота появления этого события m/N при N→∞ сходится по вероятности к среднему значению т. е. для любого > 0

lim P_r {| m/N – åP_i /N | ³ e} = 0, (4.2)

N® ¥

где m – число положительных исходов испытаний.

в)Теорема Чебышева. Если в N независимых испытаниях наблюдаются значения случайной величины X, то при N → ∞ среднее арифметическое значений сходится по вероятности к математическому ожиданию случайной величины X, т.е. для любого > 0 справедливо соотношение

lim P_r {|åx_i /N – a | ³ e} = 0, (4.3)

N® ¥

где a = M{x}, M – знак математического ожидания.

г)Обобщенная теорема Чебышева. Если – независимые случайные величины с и (дисперсии), причем все ограничены сверху одним и тем же числом, то при N → ∞ среднеарифметическое значений сходится по вероятности к среднему арифметическому их математических ожиданий, т.е. для любого имеет место

lim P_r {|åx_i /N –å a_i /N | ³ e} = 0, (4.4)

N® ¥

д)Неравенство Чебышева. Для неотрицательной функции случайной величины X и любого K > 0 имеет место неравенство

P_r{j(x) ³ K} £ M{j(x)}/K. (4.5)

В частности, если и где - среднее арифметическое, а – среднее квадратическое отклонение, то

P_r{|X - | ³ } £ 1/ . (4.6)

е)Теорема Маркова. Выражение (4.4) справедливо и для зависимых случайных величин если только

lim D{å x_i}/N ² = 0. (4.7)

N® ¥

ж) Центральная предельная теорема. Если – независимые и одинаково распределенные случайные величины с одинаковым математическим ожиданием то при N → ∞ закон распределения суммы неограниченно приближается к нормальному т.е.

lim P_r{a < (å x_i -Na)/ < b } = f(b) - f(a), (4.8)

N® ¥

где f(x) – интеграл вероятности:

f(x) = (1/

з)Теорема Лапласа. Если в каждом из N независимых испытаний события A появляется с вероятностью P, то имеет место

lim P_r{a < (m – NP)/ < b } = f(b) - f(a), (4.9)

N® ¥

где m - число реализаций события A в N испытаниях. Эта теорема является частным случаем предыдущей теоремы.

Как отмечалось выше, при статистическом моделировании основное внимание уделяется учету случайных факторов и оценке их влияния на процесс функционирования системы и ее компонентов. Для формирования на ЭВМ реализаций случайных объектов, вероятностные характеристики которых известны, сперва с помощью программных датчиков (или используя физические датчики, таблицы случайных чисел) генерируются базовые последовательности псевдослучайных чисел {r_i}, равномерно распределенных в интервале[0, 1], которые далее преобразуются в реализации других случайных объектов с заданными вероятностными характеристиками.

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

1. Моделирование случайного события. Пусть необходимо реализовать на ЭВМ случайное событие А, вероятность наступления которого равна Р_А. Стандартный алгоритм реализации этого события сводится к следующему:

- реализуем случайную величину r_i, равномерную в интервале [0, 1];

- проверяем условие r_i £ Р_А; если оно выполняется, то считаем, что событие А имеет место (т. е. оно наступило), в противном случае заключаем, что оно не наступило (вероятность этого события равна 1 - Р_А).

Основой этого простого заключения (или суждения) служит то свойство равномерного закона распределения, что вероятность попадания случайной величины в заданный интервал равна длине этого интервала, другими словами,

Р_А = . (4.10)

Итак, при выполнении неравенства r_i £ Р_А исходом испытания является событие А.

2. Моделирование полной группы событий. Техника реализации события предыдущего пункта подсказывает нам алгоритм моделирования полной группы событий А₁, А₂, … , А_n с вероятностями наступлений P₁, P₂, … , P_n, å P_k = 1. Введем обозначение

l_k = (4.11)

Тогда реализация события А_к с вероятностью Р_к будет эквивалентна выполнению условия

l_k-1 < r_i £ l_k. (4.12)

Аналогично условию (4.10), имеет место равенство

Р_k = (4.13)

Как следует из этой логики, моделирование испытания сводится к выделению интервала I_k = [l_k-1, l_k], k = 1, 2, … , n, l₀ = 0, l_п = 1, и каждому интервалу I_k будет поставлено одно определенное событие А_к, k = 1, 2, … , n.

3. Моделирование сложного события. Пусть необходимо моделировать событие С, состоящее из двух и более событий А, В, …. Для простоты рассмотрим частный случай С = (А, В), когда события А и В независимы и, пусть Р_А и Р_В – вероятности наступления этих событий соответственно. Если через и обозначить события, означающие не появление А и В соответственно, то мы будем иметь дело с моделированием возможных вариантов сложного (составного) события С: а) С = (А, В), б) С = (А, ), С = ( , В), С = ( , ), вероятности которых равны соответственно Р₁ = Р_А Р_В, Р₂ = Р_А (1 - Р_В), Р₃ = (1 - Р_А)Р_В, Р₄ = (1 - Р_А )(1 - Р_В). Эти четыре события составляют полную группу, поэтому Р₁ + Р₂ + Р₃ + Р₄ = 1.

Следовательно, для моделирования испытания можно воспользоваться приведенной выше схемой, используя для этой цели всего одну реализацию случайной величины r_i, равномерной в интервале [0, 1]:

- если 0 < r_i £ Р₁, то произошло событие (А, В);

- если Р₁ < r_i £ Р₁ + Р₂, то произошло событие (А, );

- если Р₁ + Р₂ < r_i £ Р₁ + Р₂ + Р₃, то произошло событие ( , В);

- если Р₁ + Р₂ + Р₃, < r_i £ 1, то произошло событие ( , ).

Другая модификация этого алгоритма, использующая уже две реализации случайной величины r_i и r_i+1, заключается в следующем:

а) проверить условие 0 < r_i £ Р_А; если оно выполнено, перейти к пункту б), в

противном случае – к пункту в);

б) проверить условие 0 < r_i+1 £ Р_В; если оно выполнено, то произошло

событие (А, В), в противном случае – событие (А, );

в) проверить условие 0 < r_i+1 £ Р_В; если оно выполнено, то произошло ( , В),

в противном случае – событие ( , ).

Рекомендуем студентам построить и исследовать логическую схему этого моделирующего алгоритма.

Рассмотрим теперь случай, когда события А и В статистически зависимы, и пусть заданы безусловные вероятности Р_А, Р_В и условная вероятность Р_В/А. Для моделирования нам будет необходимо также выражение для вероятности , которую можно получить из формулы для полной вероятности для события В, т. е. Р_В:

Р_В = Р_АР_В/А + (1 – Р_А) , (4.14)

откуда следует, что

= (Р_В - Р_АР_В/А)/(1 – Р_А). (4.15)

Тогда можно воспользоваться либо одной, либо двумя реализациями равномерно распределенной в интервале [0, 1] случайной величины r_i. В обоих случаях, как и выше, возможными исходами являются:

а) исход С = (А, В) с вероятностью Р₁ = Р_АР_В/А;

б) исход С = (А, ) с вероятностью Р₂ = Р_А (1 - Р_В/А);

в) исход С = ( , В) с вероятностью Р₃ = (1 – Р_А) ;

г) исход С = ( , ) с вероятностью Р₄ = (1 – Р_А)(1 - ).

При моделировании по первой схеме используется лишь одна реализация r_i для выделения (генерации) интервала, длина которого равна вероятностям Р₁, Р₂, Р₃, Р₄, причем, Р₁ + Р₂ + Р₃ +Р₄ = 1. Для второй схемы воспользуемся двумя реализациями r_i и r_i+1. Соответствующий моделирующий алгоритм имеет вид

а) проверить условие 0 < r_i £ Р_А; если оно выполнено, перейти к пункту б), в

противном случае – к пункту в);

б) проверить условие 0 < r_i+1 £ Р_В/А; ; если оно выполнено, то произошло

событие (А, В), в противном случае – событие (А, );

в) проверить условие 0 < r_i+1 £ ; если оно выполнено, то произошло

( , В), в противном случае – событие ( , ).

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

Рекомендуем студентам построить логическую схему этих моделирующих алгоритмов и изучить их внутренние связи.

4. Моделирование дискретного распределения. Пусть распределение случайной величины Х задано в виде последовательности реализаций x₁ £ x₂ £ …., £ x_n с соответствующими вероятностями p₁, p₂, … , p_n, å р_k = 1. Интегральная функция распределения этой величины Х будет определяться с помощью формулы

F_X(x) = P_r{X< x} = x_k £ x < x_k+1, k = 1, … , n. (4.16)

Этот случай эквивалентен случаю моделирования полной группы событий, рассмотренный в пункте 2, поэтому моделирующий алгоритм выглядит следующим образом. Генерируем реализацию случайной величины r_i :

- если имеет место r_i < р₁, то Х = х₁, иначе,

- если имеет место r_i < р₁ + р₂, то Х = х₂, иначе,

- ……………………………………………………………..

- если имеет место r_i < то Х = х_к,

и т. д., до значения k = n. Интересно отметить, что для этого алгоритма среднее число сравнений равно величине

M₀ = (4.17)

Рассмотрим в качестве примеров дискретного распределения следующие частные случаи, представляющие большой интерес для практики моделирования.

а) Моделирование биномиального распределения. Пусть проводятся N независимых экспериментов; в каждом эксперименте вероятность успешного исхода постоянна и равна заданной величине р. Вероятность m успешных экспериментов, m = 0, 1, … , N, определяется согласно формуле

p_m = C_N^mp^m(1 – p)^N-m, m = 0, 1, … , N, (4.18)

где C_N^m = N/m!(N – m)!. Тогда последовательность испытаний сводится к реализации следующего алгоритма:

- если имеет место r_i < р₀, то m = 0, иначе,

- если имеет место r_i < р₀ + р₁, то m = 1, иначе,

- ……………………………………………………………..

- если имеет место r_i < то m = к,

и т. д., до значения k = N. Как известно, первые два момента биномиального распределения равны M[m] = Np, D[m] = Np(1 – p), где M – знак математического ожидания, D – знак дисперсии.

б) Моделирование пуассоновского распределения. Случайная величина Х распределена по закону Пуассона, если она принимает целочисленные значения к = 0, 1, 2, … , с вероятностями

P_r{X = k} = p_k = (a^k/k!) e^-a, k = 0, 1, 2, …, (4.19)

a – среднее значение. Известно, что дисперсия этой случайной величины также равна величине a, ток что M[Х] = D[Х] = а. В принципе, моделирующий алгоритм для этого распределения можно построить по той же схеме, что и другие дискретные распределения, испытав сперва r_i < р₀, далее r_i < р₀ + р₁, и т. д. Однако существует более изящный для этого закона распределения алгоритм, который основан на двустороннем неравенстве

(4.20)

Согласно этой формуле, при нулевом начальном значении параметра k генерируются реализации случайной величины r_i, равномерно распределенной в интервале [0, 1], до тех пор, пока не будет справедливым условие (4.20). Реализующая этот алгоритм подпрограмма SUBROUTINE POISSN(a, X) на языке FORTRAN имеет вид

1. X = 0.0

2. A = exp(-a)

3. S = 1.0

4. CALL RANDUM(IX, IY, RN)

5. S = S*RN

6. IF(S – A)9,7,7

7. X = X + 1

8. GO TO 4

9. RENURN

10. END

Четвертый оператор этой подпрограммы обращается к генератору случайных чисел и формирует реализацию r_i º RN, произведение которых формируется пятым оператором. Шестой оператор проверяет знак разности S – A; если S – A < 0, управление передается оператору 10, который и возвращает текущее целочисленное значение X; если же S – A ³ 0, управление передается седьмому оператору, который прибавляет единицу к значению Х, и генерация r_i º RN повторяется.

5. Моделирование непрерывных случайных величин по методу обратной функции. Пусть необходимо генерировать на ЭВМ реализации непрерывной случайной величины Х, функция плотностираспределения которой есть f_X(x). Обозначим через F_X(x) функцию вероятности (или интегральный закон распределения), которая определяется в виде

F_X(x) = (4.21)

Сущность метода обратной функции заключается в том, что если r_i – реализация равномерно распределенной в интервале [0, 1] случайной величины R, то уравнение

r_i = F_X(x_i ) = (4.22)

определяет соответствующую реализацию x_i случайной величины Х с вероятностными характеристиками f_X(x) и F_X(x). Так как функция F_X(x) является монотонно возрастающей, ее обратная функция F_X^-1 (x) существует, поэтому из (4.22) для реализации x_i получим соотношение

x_i = F_X^-1 (r_i). (4.23)

Это преобразование и называется методом обратной функции; оно устанавливает взаимно однозначное соответствие между двумя последовательностями реализаций случайных величин {r_i} и {x_i}. Рассмотрим частные случаи.

а) Моделирование экспоненциального закона распределения. Для экспоненциального закона распределения функции f_X(x) и F_X(x) имеют вид

f_X(x) = lе^-lх, 0 £ х,

F_X(x) = 1 – е^-lх.

Решение уравнения (4.23) в этом случае равно

x_i = F_X^-1 (r_i) = -(1/l)ln(1 - r_i). (4.24)

И так как случайные числа r_i и 1 - r_i распределены одинаково, вместо (4.24), обычно используется преобразование (или моделирующий алгоритм)

x_i = -(1/l)ln r_i, i = 1, 2, … (4.25)

б) Моделирование равномерного закона распределения. Пусть необходимо сгенерировать реализации случайной величины, равномерно распределенной в интервале [a, b]. Вероятностные характеристики этой случайной величины равны соответственно

f_X(x) = 1/(b – a), a£ x£ b,

Легко проверить, что уравнению (4.23) соответствует решение

x_i = а + (b – a)r_i, i = 1, 2,… (4.26)

Согласно этой формуле, при r_i = 0 получим x_i = а, а при r_i = 1 - соответственно x_i = b.

в) Моделирование нормального закона распределения. Вероятностные характеристики нормального (или гауссова) распределения имеют вид

f_X(x) = (1/Ö2p)s_xexp(-(x – m_x)²/2s_x² ),

F_X(x) = ,

где m_x и s_x² = математическое ожидание и дисперсия случайной величины Х. Так как функция F_X(x) монотонно возрастающая, то, в принципе, для построения моделирующего алгоритма можно воспользоваться уравнением r_i = F_X(x_i ). Однако с вычислительной точки зрения более экономно построить моделирование нормально распределенной случайной величины, исходя из центральной предельной теоремы теории вероятностей, согласно которой для достаточно большого значения N величина

z = (4.27)

также распределена по нормальному закону. Практически ограничиваются значением N = 12, тогда математическое ожидание z равно 12/2 = 6, а дисперсия равна 12/12 = 1, так как величины r_i равномерно распределены в интервале [0, 1] с математическим ожиданием ½ и дисперсией 1/12.

Учитывая то обстоятельство, что величина

U = (X – m_x)/ s_x (4.28)

также распределена по нормальному закону с параметрами m_u = 0, s_u² = 1, с учетом (4.27) получим соотношение

(X – m_x)/ s_x = - 6, (4.29)

откуда получим для моделирования реализации величины Х моделирующий алгоритм в виде

x = m_x + s_x ( - 6). (4.30)

г) Моделирование закона распределения с кусочно- постоянной аппроксимацией.

Когда случайная величина имеет довольно сложный закон распределения, удобно аппроксимировать функцию плотности распределения с помощью кусочно-постоянной функцией, что, в принципе, позволяет получить универсальный моделирующий алгоритм. Пусть функция распределения f_X(x) определена на интервале [a, b], который разбит на l подинтервалов, внутри которых функция f_X(x) постоянна и равна c_k, k = 1, … , l. Именно в такой форме представляются результаты многочисленных экспериментов, в том числе и экспериментов с моделями систем, так что подобные алгоритмы достаточно универсальны и широко применяются в практике исследования и моделирования. Пусть далее r_к^* - реализации случайной величины, равномерной в каждом из этих подинтервалов. Тогда величина x = a_k-1 + r_к^* будет распределена по закону f_X(x), где a_k-1 – левый край k - го подинтервала. Алгоритм реализации случайной величины r_к^*, равномерной в подинтервале [a_k-1, a_k], был рассмотрен нами выше.

Для практических целей более удобно аппроксимировать закон распределения f_X(x) таким образом, чтобы вероятность попадания в любой участок [a_k-1, a_k] была постоянна и не зависела от номера этого участка. Такое предположение означает, что границы этих участков удовлетворяют уравнению

(4.31)

причем a₀ = a, a_l = b. Соответствующий моделирующий алгоритм получается в виде следующей последовательности

1. получить реализацию r_i, равномерную в интервале [0, 1];

2. с помощью этой реализации r_i случайным образом выбрать подинтервал [a_k-1, a_k];

3. генерируется реализация r_i+1 и масштабируется с целью приведения ее к подинтервалу [a_k-1, a_k] путем преобразования r_i+1: = (a_k - a_k-1) r_i+1;

4. вычисляется реализация x_i = a_k + r_i+1.

Последовательность {x_i} оказывается распределенной по закону f_X(x).

д) Моделирование векторной случайной величины Z = (X, Y). В предположении, что закон распределения f_Z(x, y) задан, реализации векторной величины z можно построить на основе следующих соображений. Вычислим закон распределение f_X(x) путем интегрирования

f_X(x) = (4.32)

и определим функцию условного распределения f_Z(y/х), согласно выражению

f_Z(y/х) = f_Z(x, y)/ f_X(x). (4.33)

Тогда вектор Z = (X, Y) можно смоделировать следующим образом:

1. получить реализацию r_i, равномерную в интервале [0, 1], и с помощью функции F_X (x) смоделировать реализацию х = F_X^-1 (r_i);

2. получить реализацию r_i+1, равномерную в интервале [0, 1], и с помощью функции F_Z(y/х) смоделировать реализацию y; тогда пара реализаций (х, y) = z будет распределена по закону f_Z(x, y) = f_Z(y/х)f_X(x).

Функции F_X(x) и F_Z(y/x) находятся из интегральных выражений

F_X(x) = F_Z(y/x) = . (4.34)

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

6. Моделирование марковского процесса. В разделе 1.2 уже шла речь о дискретных цепях Маркова, которые составляют важный класс дискретных управляемых вероятностных процессов. Напомним, что в любой дискретный момент времени (или в любом такте) автомат (или управляемая система) может находиться в одном из состояний z_k, k = 1, … , K. Матрица переходов из одного состояния в другое состояние задается в виде

, 0 £ p_ij £ 1, i, j = 1, … , K. (4.35)

Вероятности p_ij интерпретируются как переходные вероятности из состояния z_i , в котором находится система (процесс, автомат и т. д.) в момент времени t, в состояние z_j в следующий момент времени t + 1. Предполагается, что переходные вероятности удовлетворяют условию нормировки

(4.36)

Если задаваться также начальными вероятностями p_i(0), i = 1,2, … , K, и обозначить через p_i(п), i = 1,2, … , K, вероятности того, что после п переходов система будет находиться в состоянии z_i, тогда моделирование процесса можно осуществить на основе так называемого событийного подхода, который каждому состоянию z_i ставит в соответствие событие А_i, i = 1, … , K.

Первое испытание проводим с помощью реализации r_i, используя при этом начальные вероятности p_i(0), i = 1,2, … , K, сумма которых равна единице. Этим самым определяется событие, обозначим его через , которое является исходом первого испытания. Далее с помощью реализации r_i+1 проводим второе испытание, используя вероятности , j = 1, … , K, первый индекс которых соответствует номеру события , и обозначим исход через, A_i2, и т. д. Продолжая эти испытания за п шагов, получим искомый процесс. Иногда в качестве начальных вероятностей выбираются p_i(0) = 1/К, i = 1,2, … , K.

7. Псевдослучайные числа и алгоритмы их генерации. Выше мы отметили, что при применении метода статистического моделирования используются программные датчики для получения последовательности реализации случайной величины, равномерно распределенной в интервале [0, 1]. Функция плотности распределения такой величины постоянна в этом интервале и также равна единице, а функция вероятности (интегральный закон распределения) определяется соотношением

(4.37)

Ранее отмечалось также, что математическое ожидание этой величины равно ½, а дисперсия равна 1/12. Получить на ЭВМ реализации непрерывной случайной величины не представляется возможным из-за дискретности самого устройства, возможно лишь получить дискретные значения каждой реализации, которые в случае k – разрядной машины имеют вид

r_i = i/(2^k – 1), i = 0, 1, 2, …2^k – 1. (4.38)

Если все эти значения имеют одинаковую вероятность, равную

p_i = 1/2^k, i = 0, 1, 2, …2^k – 1, (4.39)

тогда первые два момента этой дискретной случайной величины равны соответственно

M[R] = (4.40)

D[R] =

= (4.41)

Из этих соотношений видно, математические ожидания непрерывной и дискретной величин совпадают, однако их дисперсии отличаются друг от друга множителем

u(k) = (2^k + 1)/(2^k – 1), (4.42)

который стремится к единице при достаточно большом значении k.

Программная реализация дискретных случайных величин, равномерно распределенных в интервале [0, 1] (так называемых псевдослучайных чисел, ПСЧ) основана на известном из теории чисел свойстве конгруэнтности: числа a и b называются конгруэнтными (сравнимыми) по модулю m, или в виде формулы

a º b (mod m), (4.43)

если их разность кратна числу m, т. е. можно записать a - b = km, где k – некоторое целое число.

В технике моделирования широкое распространение получили следующие алгоритмы, действие которых основано на конгруэнтности чисел:

а) мультипликативный конгруэнтный алгоритм, основанный на формуле

X_i+1 º lX_i (mod M), (4.44)

б) смешанный конгруэнтный алгоритм, основанный на формуле

X_i+1 º lX_i + m (mod M), (4.49)

в) аддитивный конгруэнтный алгоритм, основанный на формуле

X_i+1 º (X_i + X_i-1) (mod M). (4.50)

Поясним смысл алгоритма (4.44). Числа X_i , l и M являются неотрицательными целыми числами.В двоичном счислении умножаются числа X_i и l, при этом полученное новое целое число имеет 2b разрядов, где b – длина машинного слова(число двоичных цифр (бит) в машинном слове); далее b старших разрядов автоматически отбрасываются, и оставшиеся b младших разрядов принимается за следующее целое число X_i+1, и т. д. Если эти числа разделить на самое большое число, которое можно получить на компьютере (в данном случае – это число 2^b), получим случайные реализации r_i, равномерные в интервале [0, 1]. Действие остальных алгоритмов аналогично этому алгоритму.

При успешном выборе начального значения X₀, а также величин l и M достигается приемлемое качество сгенерированной с помощью конгруэнтных алгоритмов последовательности псевдослучайных чисел {r_i}. Так, например, для двоичной системы счисления величина M выбирается равной 2^b, величина l, в свою очередь, выбирается из правила l = 8Т ± V₀, где Т – любое целое положительное число, V₀ – любое положительное, но нечетное число (например, V₀ = 3) . При таком выборе максимальный период последовательности псевдослучайных чисел {r_i} получается равным

P = 2^b-2 = M/4, b > 2. (4.51)

Большую длину последовательности {r_i} можно получить, частично пожертвовав скоростью вычислений: оказывается, что если выбрать величину M, равную наибольшему простому числу, которое меньше чем 2^b, а величину l - равной первообразному корню из M, максимальная длина последовательности {r_i} будет увеличена от M/4 до значения M – 1. Например, при b = 36 наибольшее простое число, меньшее 2³⁶, есть 2³⁵ – 31, следовательно, M = 2³⁵ – 31 = 34 359 738 337. Первообразные корни из M равны 5⁵ и 5¹⁵, поэтому можно брать l = 5⁵, что, впрочем, также согласуется с формулой l = 8Т ± V₀ при V₀ = 3.

Так как эффективность метода статистического моделирования и достоверность полученных результатов в значительной мере зависят от качества последовательности {r_i}, ее подвергают тщательной проверке. В практике статистического моделирования систем наиболее часто используемыми являются следующие методы: а) проверка (или тест) на равномерность; б) проверка на стохастичность; в) проверка на независимость. Рассмотрим вкратце основную сущность этих процедур.

а) Проверка на равномерность. Суть этой проверки сводится к тому, чтобы установить статистическую значимость гипотезы о том, что эмпирическая кривая распределения (гистограмма) близка к равномерному закону распределения на интервале [0, 1]. Пусть этот интервал разбит на l равных частей, тогда при генерации последовательности {r_i} каждая реализация попадает в один из подинтервалов единичного отрезка с одинаковой вероятностью, равной p_k = 1/l, k = 1, … , l. Обозначим через m_i, i = 1, … , l, количество попаданий в каждый из этих подинтервалов, так что å m_i = N, где N = общее количество реализации в последовательности {r_i}. Величины m_i/N, i = 1, … , l, характеризуют относительные частоты попадания в каждый из подинтервалов. Поэтому, в принципе, справедливость гипотезы о близости теоретического и эмпирического распределений можно установить на основе критерия согласия «Хи – квадрат» Пирсона

(4.52)

Часто для этой цели используется косвенный метод проверки, суть которого заключается в следующем. Последовательность {r_i} разбивается на две подпоследовательности {r_1, r₃, … , r_2i-1} и {r₂, r₄, … , r_2i}, i = 1, … , N. Если после N/2 испытаний, когда сгенерировано N реализаций, условие

r²_2i-1 + r²_2i £ 1, i = 1, …, N (4.53)

выполняется ровно k раз, k £ N/2, то, согласно предельной теореме теории вероятностей, отношение k/N/2 при достаточно большом значении N должностремиться к величине p/4, поскольку, по геометрическим соображениям, точка (r_2i-1, r_2i) всегда принадлежит единичному квадрату с площадью S_кв = 1, а условие (4.49) устанавливает факт попадания этой же точки в четверть единичного круга с площадью S_кр = pr²/4 = p/4, так что отношение площади четверти единичного круга к площади единичного квадрата равно S_кр / S_кв = p/4.

б) Проверка на стохастичность. Эту проверку можно осуществить, используя так называемые комбинации и серии частных подпоследовательностей реализации псевдослучайных чисел {r_i}. Суть проверки на основе комбинаций сводится к определению закона распределения длин участков между единицами (или нулями) в k – разрядных двоичных числах X_i, i = 1, 2, … , N, которые получаются алгоритмически на основе формул (7.8) – (7.10), где N – длинапоследовательности {r_i}. При достаточно большом значении N обычно проверяются все k разрядов или только l старших разрядов случайных чисел X_i. Так как вероятность появления единицы или нуля в каждом разряде этих чисел одинакова и равна p = 0.5, то вероятность появления m единиц в l разрядах двоичных чисел X_i будет определяться на основе формулы для биномиального распределения

p(m, l) = C_l^mp^m(l)(1 – p(1))^l-m = C_l^mp^l(l), (4.54)

где p(1) = p(0) = 0.5, C_l^m = l!/m!(l – m)!. Тогда в последовательности длины N ожидаемое число появления чисел X_i с m единицами в проверяемых l разрядах будет равно величине

N_m = N p(m, l) = N C_l^mp^l(l). (4.55)

Имея теоретические и эмпирические значения вероятностей p(m, l) или чисел N_m, можно проверить стохастичность последовательности {r_i} с помощью некоторого критерия согласия.

Проверка на основе использования серий сводится к тому, что последовательность {r_i} разбивается на две группы таким образом, чтобы выполнялось условие

x_i = (4.56)

В этом выражении величина p из интервала (0, 1). Под серией подразумевается любой отрезок последовательности {r_i}, состоящий из следующих друг за другом элементов типа a или b, например, ряд a a b b b a a a a a a b b a a a a a b b b b aaa… состоит из семи серий, длина которых равна соответственно числам 2, 3, 6, 2, 5, 4, 3. Очевидно, что вероятность появления серии длиной m в последовательности длиной l в N реализациях также будет определяться на основе формулы (7.14), причем p(1) = p(0) = p, m = 0, 1, … , l; l = 1, 2, … , k. Сравнивая эмпирические и теоретические вероятности на основе некоторого критерия согласия, как и выше, можно судить о стохастичности последовательности {r_i}.

в) Проверка независимости элементов последовательности. Независимость элементов последовательности {r_i} можно установить с помощью корреляционных отношений. Так, например, образовав на основе последовательности {r_i} другую последовательность {v_i} на основе соотношения v_i = r_i+k, где k – величина сдвига последовательности, и вычислив корреляционный момент K_rv для этих двух рядов и нормализовав последний путем разделения на произведение среднеквадратических отклонений и , т. е. вычислив величину r_rv = K_rv/ , независимость элементов последовательности можно проверить с помощью неравенства

r_rv ² £ b²/N, (4.57)

где b - вероятность доверия, N – длина ряда {r_i}.

При выполнении расчетов на ЭВМ в качестве оценок ² и ² используются соответствующие эмпирические дисперсии

D(r_i) = (1/(N – k)) , (4.58)

D(r_i+k) = (1/(N–k)) , (4.59)

тогда эмпирический корреляционный момент оценивается по формуле

r_rv (k) = ((1/(N – k)) .(4.60)

При генерировании последовательностей {r_i} и {v_i} формирование величин år_i, år_i+k, år_ir_i+k, år_i² и år_i+k² происходит автоматически в соответствующих ячейках памяти.

Важными характеристиками последовательности {r_i} являются длина периода Р и длина отрезка апериодичности L. Длина отрезка апериодичности L есть наибольшее целое число, такое, что при 0 £ j < i £ L событие A: = {r_i = r_j} не имеет места, другими словами, все числа r_i в пределах отрезка апериодичности не повторяются. Величины Р и L устанавливаются экспериментально. Очевидно, что если при статистическом моделировании длина последовательности {r_i} превышает величину L, то испытания могут повторяться при тех же условиях, и при этом новых статистических результатов не будет получено.

Отметим, наконец, что все названные выше критерии оценки качества последовательности {r_i} являются лишь необходимыми условиями; их достаточность в каждом конкретном случае подлежит дополнительному исследованию.