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

1.6. Метод статистического моделирования На этапе исследования и проектирования систем при построении и реализации машинных моделей широко используется метод статистических испытаний (метод Монте-Карло), который базируется на использовании случайных чисел. Метод статистических испытаний – это метод решения невероятностной проблемы вероятностным способом. Часто его называют методом Монте-Карло, по названию проекта США по ядерной технике, где он был впервые предложен. Явное представление времени здесь отсутствует. Суть метода в том, что процесс описывают формулами и логическими выражениями на ЭВМ. Затем в модель вводят случайно изменяющиеся факторы и оценивают их влияние на показатели процесса. Результаты оценки подвергают статической обработке. Идею этого метода покажем на примере расчета площади криволинейной фигуры (рис. 2). Рис. 2. Расчет площади криволинейной фигуры:а – аналитическая модель; б – метод Монте-Карло При традиционном подходе строят аналитическую зависимость Y=φ(X), по которой рассчитывают площадь фигуры (рис. 2). На этом прямоугольнике в случайном порядке разбрасывают точки. Затем оценивают, какая доля всех точек попала внутрь фигуры. Эта доля соответствует доле площади фигуры от рассчитанной площади прямоугольника. Чем больше точек, тем точнее можно определить площадь фигуры. Метод статистических испытаний, как первый и широко распространившийся метод имитационного моделирования, часто вообще считают единственным методом имитации систем. Динамику системы описывают на универсальном языке программирования в виде последовательности уравнений с детерминированными –Х и случайными ~X(Δt) коэффициентами (рис. 3). Время моделирования разбивают на одинаковые шаги Δt. EMBED Word.Picture.8 Рис. 3. Имитация поведения системы методом статистических испытаний На каждом шаге Δt изменяют значения случайных коэффициентов ~X(Δt), для которых по уравнениям рассчитывают изменения выходной величины ~Y(Δt). Каждый эксперимент представляет собой расчет уравнений с шагом Δt. В результате выявляют связь выходных величин с входными величинами. Поскольку шаг времени Δt выбирают для самого быстрого элемента системы, то в результате моделирования получают много ненужной информации для шагов, когда никаких изменений этого элемента в системе не происходит. Кроме того, разработка имитационной модели на универсальном языке программирования требует специальной квалификации и часто остается без изменений в системе. Существуют две области применения метода статистического моделирования: для изучения стохастических систем; для решения детерминированных задач. В настоящее время моделирование по методу Монте-Карло широко применя­ется при решении определенных задач статистики, которые не поддаются анали­тической обработке. Этот тип моделирования применялся для оценки критиче­ских значений или достоверности критерия при проверке гипотезы.^ 1.6.1. Общая структура статистической модели Задачи статистического моделирования: 1) построение объекта моделирования; 2) формирование случайных взаимодействий; 3) организация статистической обработки данных моделирования; 4) задача планирования эксперимента. На рис. 4 приведена общая структура статистической модели. Блоки I, II, III соответствуют первым трем задачам. Рис. 4. Общая структура статистической модели^ 3. Моделирование случайных процессов Еще раз отметим, что имитационная модель позволяет исследовать поведе­ние различных систем с учетом влияния случайных факторов. Эти факторы в зависимости от их природы могут быть отражены в модели как случайные события, случайные величины (дискретные или непрерывные) или как случай­ные функции (процессы). Например, если с помощью создаваемой имитационной модели предполагается исследовать надежность вычислительной системы, то возникновение отказа будет представлено в модели как случайное событие. Если же модель предназначена для оценки временных параметров процесса обслуживания клиентов в автомастерской, то интервал времени до появления очередного клиента удобнее всего описать как случайную величину, распределённую по некоторому закону. К общим принципам имитации случайных воздействий на ЭВМ относятся: 1) Формирование базовой случайной величины (CB). 2) Преобразование базовой случайной величины в значения случайных величин, распределенных по требуемому закону.^ 3.1. Способы формирования базовой случайной величины В основе базовой СВ обычно используется СВ равномерно распределенная на интервале [0,1]. Рассмотрим общий случай распределения СВ на интервале [a,b]. Непрерывная СВ имеет равномерное распределение в интервале [a,b], если ее функции плотности (рис. 13) и распределения (рис. 14) соответственно примут вид: Рис. 13. Функция плотности для равномерного распределения Рис. 14. Функция равномерного распределения Определим числовые характеристики случайной величины , принимающей значения x: математическое ожидание (уравнение 17), дисперсию (уравнение 18) и среднее квадратическое отклонение (уравнение 19):, (17), (18). (19) При моделировании систем на ЭВМ приходится иметь дело со случайными числами интервала [0,1], когда границы интервала a=0, b=1. Рассмотрим частный случай равномерного распределения, когда функции плотности и распределения соответственно имеют вид: ,.Математическое ожидание такого распределения – , дис­персия – и среднеквадратическое отклонение – . Это распределение требуется получить на ЭВМ. Но полу­чить его на цифровой ЭВМ невозможно, т.к. машина оперирует с n-раз­рядными числами. Поэтому на ЭВМ вместо непрерывной совокупности равно­мерных случайных чисел интервала [0,1] используют дискретную по­следовательность 2n случайных чисел того же интервала. Закон рас­пределения такой дискретной последовательности называют квазиравномерным распределением. ^ 3.2. Получение квазиравномерных чисел Случайная величина , имеющая квазиравномерное распределение в интервале [0,1], принимает значения с вероятностями . Математическое ожидание и дисперсия квазиравномерной СВ соответственно имеют вид:,. В первом случае используем соотношение:. (20) Во втором случае имеем соотношение:. (21) Таким образом, математическое ожидание квазиравномерной случайной величины совпадает с математическим ожиданием равномерной случайной последовательности интервала [0,1], а дисперсия отличается множителем (2n + 1) / (2n – 1), который при достаточно больших n близок к единице. На ЭВМ невозможно получить идеальную последовательность случайных чисел, т.к. можно оперировать только с конечным множеством чисел, и для получения значений х случайной величины используют формулы; поэтому такие последовательности называют псевдослучайными.^ 3.3. Моделирование случайных событий Для моделирования случайного события А, вероятность которого равна Рс, достаточно сформировать одно число r, равномерно распределенное на интервале [0,1]. При попадании r в интервал [0, Рс] считают, что событие А наступило, в противном случае – не наступило, т.е.:. (22) Пусть, например, вероятность отказа вычислительной системы составляет 0,3. Чтобы определить, возникнет ли он на очередном шаге моделирования, достаточно сгенерировать с помощью датчика одно случайное число r и срав­нить его с вероятностью отказа (рис. 15). а) б) Рис. 15. Моделирование случайного события: а) отказ «произошёл»; б) отказ «не произошёл» Для моделирования полной группы N несовместных событий А={A1, A2, …,AN} с вероятностями соответственно P1, Р2,..., PN также достаточно одного значения r: событие Аi из группы А считается наступившим, если выполняется условие:. (23) Предположим, что в каждый момент времени может происходить обращение только к одному из трех модулей оперативной памяти вычислительной системы. Вероятности обращения к каждому из них Р1, Р2 и Р3 равны соответственно 0.3, 0.5 и 0.2. Чтобы узнать, из какого именно модуля будут считаны данные, необходимо определить, в какой интервал попадает полученное от датчика случайное число r (рис. 16). EMBED Word.Picture.8 Рис. 16. Моделирование полной группы из трех несовместных событий Если группа событий А не полна, то вводят фиктивное событие AN+1 с вероятностью PN+1 такой, что сумма вероятностей становится равной 1:, т.е. дополняют группу ^ А до полной. (24)После этого генерируют число r и проверяют условие (23). При А = AN+1считают, что ни одно событие из исходной группы А не наступило. Для имитации зависимых событий А и В (В зависит от А) необходимо знать безусловную вероятность Р(А) события А и условные вероятности Р(В/А) и Р(В/А). Сначала описанным выше способом имитируется появление собы­тия А, в зависимости от исхода выбирается одна из вероятностей Р(В/А) или Р(В/А), и по той же технологии определяется наступление события В.^ 3.4. Моделирование случайных величин При моделировании дискретных случайных величин наиболее часто ис­пользуются два метода: метод последовательных сравнений; метод интерпретации.^ 3.4.1. Метод последовательных сравнений Число r последовательно сравнивают со значением суммы Р1+Р2+..., где Р1 – вероятность наименьшего значения СВ Y, Р2 – вероятность второго по величине значения. При первом выполнении условия r> проверка прекращается и дискретная СВ Y считается принявшей значение уi. Процесс можно ускорить, применяя методы оптимизации перебора: дихо­томии, ранжирования Р и т.д. Величины Pi рассчитывают по функциям распределения вероятности, соот­ветствующим моделируемому закону.^ 3.4.2. Метод интерпретации Метод основан на физической трактовке моделируемого закона распреде­ления. Например, биномиальное распределение описывает число успехов в n независимых испытаниях с вероятностью успеха в каждом испытании Р и ве­роятностью неудачи g = 1 – Р. При моделировании этого распределения с помощью метода интерпрета­ции выбирают n независимых СЧ, равномерно распределенных на интервале [0, 1], и подсчитывают количество тех из них, которые меньше Р. Это число и является моделируемой СВ.^ 3.5. Моделирование непрерывных случайных величин При моделировании непрерывных СВ с заданным законом распределения могут использоваться три метода: нелинейных преобразований; композиций; табличный. Первые два метода требуют от разработчика модели весьма серьезной математической подготовки и в данном пособии не рассматриваются. Третий метод, условно названный «табличным», основан на замене закона распределения не­прерывной СВ специальным расчетным соотношением, которое позволяет вы­числять значение СВ по значению случайного числа, равномерно распределенного все на том же интервале [0, 1]. Такие соотношения получены практически для всех наиболее распространенных видов распределений и приведены в справочной литературе. В качестве примера ниже приведено расчетное соотноше­ние для показательного закона распределения показательного:, (25) где – параметр показательного распределения,r – равномерно распределённое СЧ. Технологию моделирования непрерывных случайных величин поясним на примере. Специалисты в области надежности технических систем доказали, что интервалы времени между отказами распределены по экспоненциальному зако­ну. Чтобы получить в модели величину промежутка времени между двумя со­седними отказами, достаточно сгенерировать случайное число r и подставить его в выражение (25); разумеется, при этом должен быть известен параметр распределения . В одной и той же имитационной модели могут фигурировать различные случайные факторы, одни могут быть представлены как случайные события, другие – как случайные величины; более того, если моделируется поведение достаточно сложной системы, то ее функционирование может быть связано с возникновением нескольких типов событий и учетом большого числа случай­ных величин, распределенных по различным законам. Если же моделирование всех случайных факторов основано на использовании одного датчика, генерирующего одну «общую» последовательность случайных чисел, то с математической точки зрения их нельзя считать независимыми. В связи с этим для моделирования каждого случайного фактора стараются использовать отдель­ный генератор, или, по крайней мере, обеспечивать создание новой последова­тельности случайных чисел. Во многих специализированных языках и пакетах моделирования такая возможность предусмотрена.^ 3.6. Способы получения случайных чисел Моделирование любой системы или процесса, содержащих случайные компонен­ты, предполагает использование метода генерирования чисел, которые в опреде­ленном смысле являются случайными. Существуют два способа формирования случайных чисел (CЧ): 1) с помощью специального физического датчика; 2) с помощью специальных программ.^ 3.6.1. Физический способ Для более быстрого генерирования слу­чайных чисел были созданы механизированные устройства: в конце 1930-х годов Кендалл и Бабингтон-Смит использовали быстро вращающийся диск для подготовки таблицы, содержащей 100000 случайных однозначных чисел. Позднее разработаны электрические схемы, основанные на произвольно пульсирующих электрон­ных лампах, которые выдавали до 50 случайных чисел в секунду. Одна из таких схем реализована в электронном индикаторном устройстве для получения слу­чайных чисел, исполь­зовавшемся Британским почтовым ведомством для выбора победителей в лотерее. Другое электронное устройство применялось компанией Rand Corporation для генерирования таблицы, содержащей 1000000 случайных цифр. Было изобретено и множество других схем, таких как произвольный выбор чисел из телефонных книг, результатов переписи населения или использование цифр в числе до 100000 десятичных разрядов [8]. С распространением компьютеров (и моделирования) все более пристальное внимание стало уделяться методам генерирования, или генераторам случайных чисел, совместимых со способом работы компьютеров. При этом способе генерации случайные числа вырабатываются специальной электронной приставкой – генератором (датчиком) СЧ, служащей в качестве одного из внешних устройств ЭВМ. В качестве физического эффекта, лежащего в основе таких генераторов чисел, чаще всего используются шумы в электронных и полупроводниковых приборах, явления распада радиоактивных элементов и т.д. (количество радиоактивных частиц). Преимуществом физических датчиков является высокая скорость формирования случайных чисел. К недостаткам физических датчиков случайных чисел относятся: 1) изготовление отдельного прибора; 2) не позволяют гарантировать качество последовательности непосредственно во время моделирования системы на ЭВМ, а также повторно получать при моделировании одинаковые последовательности чисел.^ 3.6.2. Программный способ Наибольшее применение в практике моделирования на ЭВМ для генерации последовательностей псевдослучайных чисел получили алгоритмы вида: xi+1 = Ф(xi), представляющие собой рекуррентные соотношения первого порядка, для которых начальное число x0 и постоянные параметры заданы. Хороший арифметический генератор случайных чисел должен обладать следу­ющими свойствами. 1. Получаемые числа должны быть равномерно распределены в интервале [0,1] и не должны иметь корреляции друг с другом, в противном случае результаты моделирования могут оказаться полностью недействительными. 2. Чтобы генератор можно было использовать на практике, он должен обладать быстродействием и не требовать больших затрат памяти. 3. Генератор должен обеспечивать возможность точно воспроизводить заданный поток случайных чисел. 4. В генераторе должен быть предусмотрен простой способ получения отдель­ных потоков случайных чисел. Поток – это про­сто часть последовательности случайных чисел, производимых генератором, очередной поток начинается в том месте, где заканчивается предыдущий.^ 3.6.2.1. Метод серединных квадратов Одной из исторически первых процедур получения псевдослучайных чисел была процедура, получившая название серединных квадратов (была предложена фон Нейманом и Метрополисом в 1940-х годах). Алгоритм получения последовательности случайных чисел методом серединных квадратов сводится к следующему: Пусть имеется 2n-разрядное число, меньше 1:. Возведем его в квадрат:, а затем возьмем средние 2n разрядов:, которые и будут очередным числом.Например: и т.д. Недостатком этого метода является наличие корреляции между числами последовательности, а в ряде случаев случайность вообще может отсутствовать.Например: Метод срединных квадратов вовсе не является случайным, то есть непредсказуемым (это наиболее существенный его недостаток). На самом деле, если мы знаем одно число, следующее число является полностью предопределенным, по­скольку правило его получения неизменно. Фактически, когда задается х0, предопределяется вся последовательность чисел xi. Это замечание касается всех арифметических генераторов.^ 3.6.2.2. Линейные конгруэнтные генераторы Широкое применение при моделировании на ЭВМ получили линейные конгруэнтные процедуры генерации псевдослучайных последовательностей. В них последовательность целых чисел х1, х2,… определяется по рекурсивной формуле, (26) где m (модуль), λ (множитель), c (приращение) и х0 (начальное число или значение) являются неотрицательными целыми числами. Таким образом, согласно формуле (26) для получения xi нужно разделить axi-1 + c на m, то есть xi будет ос­татком этого деления. Два целых числа α и β конгруэнтны (сравнимы) по модулю m (m – целое число), если α – β делится на m без остатка и если числа α и β дают одинаковые остатки от деления на m.Например: 125 и 5 125 ≡ 5 (mod 10) m = 10. Конгруэнтная процедура получения последовательностей случайных квазиравномерно распределенных чисел может быть реализована мультипликативным (с параметром c= 0) либо смешанным методом (с параметром c> 0).^ 3.6.2.3. Мультипликативный метод Задает последовательность неотрицательных целых чисел {xi}, не превосходящих m по формуле:. (27) Для машинной реализации m = pq, где p = 2, а q – число бит в машинном слове (286 – 16; 386, 486 – 32; Pentium – 64).Алгоритм построения последовательности сводится к следующим шагам: 1) выбрать в качестве x0 произвольное нечетное число; 2) вычислить коэффициент λ = 8t ± 3, где t – любое целое положительное число; 3) найти произведение λ^ X0, содержащие не более 2q значащих разрядов; 4) взять q младших разрядов в качестве первого члена последовательности X1; 5) определить дробь x1 = X1/2q – это и есть искомое число; 6) присвоить X0 = X1; 7) вернуться к 3 пункту.^ Пример: Получим последовательность чисел для q = 4 1) X0 = 7; 2) t = 1, λ = 11 или 5 (берем 5); 3) λX0  5  7 = 35  001000112; 4) берем 4 младших разряда 00112(310); 5) x1 = 3/16 = 0,1875 – первое число; 6) ^ X0 = 00112; Для второго числа X0 = 112 (3) 3) λX0  5  3 = 15  000011112; 4) X1 = 11112 (1510); 5) x1 = 15/16 = 0,9375 – второе число; 6) ^ X0 = 11112; Для третьего числа X0 = 11112; 3) λX0  5  15 = 75  010010112; 4) X1 = 10112 (1110); 5) x1 = 11/16 = 0,6875 – третье число; 6) X0 = 1011 и т.д. ^ 3.6.2.4. Cмешанный метод Сначала выбирается значение m. Чтобы получить длинный период и высокую плотность величин xi в интервале [0, 1], величина m должна иметь большое значение. Самым удачным выбором является m=2b, где b – число битов в слове задействованного компьютера.^ 3.7. Проверка качества последовательностей псевдослучайных чисел Эффективность статистического моделирования систем на ЭВМ и достоверность получаемых результатов зависит от качества исходных (базовых) последовательностей псевдослучайных чисел, которые являются основой для получения стохастических воздействий на элементы моделируемой системы. Поэтому перед моделированием на ЭВМ необходимо убедиться в том, что исходная последовательность псевдослучайных чисел удовлетворяет предъявляемым к ней требованиям. При моделировании важными характеристиками качества генератора являются длина периода p и длина отрезка апериодичности L. Длина отрезка апериодичности L псевдослучайной последовательности есть наибольшее целое число, такое, что все числа xi в пределах этого отрезка не повторяются. Способ экспериментального определения длины периода p и длины отрезка апериодичности L сводится к следующему: 1. Запускается программа генерации последовательности {xi} с начальным значением x0 и генерируется V чисел xi. В большинстве случаев V – (1 – 5)106. Генерируются числа и фиксируется число xV. 2. Затем программа запускается повторно с начальным числом x0 и при генерации очередного числа проверяется истинность события p{xi = xV}. Если это событие истинно i = i1 и i = i2 (i1 3. Проводится запуск программы генерации с начальными числами x0 и xp. При этом фиксируется минимальный номер i = i3, при котором истинно событие , и вычисляется длина отрезка апериодичности L = i3 + p. Рис. 17. Представление определения длин периода P и отрезка периодичности L Применяемые в имитационном моделировании генераторы случайных чисел должны пройти тесты на пригодность. Существуют тесты двух типов. ^ Эмпирические тесты – это обычный тип статис­тических тестов, они основаны на действительных значениях xi, выдаваемых генератором. Теоретические тесты не являются тестами в том смысле, в каком они предусматриваются в статистике. Однако в них используются числовые парамет­ры, чтобы оценить генератор глобально без фактического генерирования некото­рых или всех значений xi. Основные анализируемые характеристики генерируемых датчиком последовательностей: равномерность; стохастичность (случайность); независимость. Рассмотрим методы проведения такого анализа, наиболее часто применяемые на практике.^ 3.7.1. Проверка равномерности Проверка равномерности может быть выполнена с помощью гистограммы относительных частот генерируемой случайной величины. Для ее построения интервал (0;1) разбивается на m равных частей и подсчитывается относительное число попаданий значений случайной величины в каждый интервал. Чем ближе огибающая гистограммы к прямой, тем в большей степени генерируемая последовательность отвечает требованию равномерности распределения (рис. 18). EMBED Word.Picture.8 Рис. 18. Гистограмма относительных частот случайных чисел^ 3.7.2. Проверка стохастичности Рассмотрим один из основных методов проверки – метод комбинаций. Суть его сводится к следующему. Выбирают достаточно большую последовательность случайных чисел xi и для нее определяют вероятность появления в каждом из xi ровно j единиц. При этом могут анализироваться как все разряды числа, так и только l старших. Теоретически закон появления j единиц в l разрядах двоичного числа может быть описан как биномиальный закон распределения (исходя из независимости отдельных разрядов). Тогда при длине выборки ^ N ожидаемое число появлений случайных чисел xi с j единицами в проверяемых l разрядах будет равно:, (28) где – число комбинаций (сочетаний) j единиц в l разрядах,рl (1) – вероятность появления единицы в двоичном разряде, р(1)=0.5. Для полученной последовательности определяется эта же характеристика. Проверка соответствия реального значения теоретическому выполняется с по­мощью одного из статистических критериев согласия.^ 3.7.3. Проверка независимости Проверка независимости проводится на основе вычисления корреляцион­ного момента. Напомним, что две случайные величины а и b называются независимыми, если закон распределения каждой из них не зависит от того, какое значение приняла другая. Для независимых случайных величин корреляционный момент равен нулю. Для оценки независимости элементов последовательности поступают сле­дующим образом. Вводят в рассмотрение дополнительную последовательность Y, в которой yi = xi + t, где t – величина сдвига последовательности Y относительно исход­ной последовательности X.Вычисляют коэффициент корреляции случайных величин X и Y, для чего используются специальные расчетные соотношения.