маткад. МИМЭП_лаб_практикум в Mathcad. Математическое и имитационное моделирование экономических процессов
 Скачать 4.96 Mb.
|
|
7. Лабораторная работа №7. Генераторы случайных величин с равномерным распределением 7.1. Общие сведения Для построения имитационных моделей необходимо иметь возможность генерирования случайных величин либо с помощью таблиц, либо по теоретическим законам распределения вероятностей с требуемыми параметрами. Для этой цели используются случайные числа или выборки по методу Монте- Карло. Если имитационная модель просчитывается на ЭВМ, мы должны иметь возможность: 1) получать равномерно распределенные случайные числа в интервале [0; 1]; 2) использовать эти случайные числа для генерации случайных величин с требуемыми характеристиками. Библиотеки программ всех ЭВМ включают с этой целью специальные стандартные подпрограммы. С помощью рекуррентных математических методов реализовано несколько алгоритмов генерирования псевдослучайных чисел. Мы называем эти числа псевдослучайными потому, что фактически они, даже пройдя все статистические тесты на случайность и равномерность распределения, остаются полностью детерминированными. Это значит, что если каждый цикл работы генератора начинается с одними и теми же исходными данными (константами и начальными значениями), то на выходе мы получаем одинаковые последовательности чисел. 7.2. Моделирование случайных величин с равномерным распределением в интервале [0; 1] Плотность вероятности ( ) f x случайной величины x , равномерно распределенной в интервале [0; 1], равна 1, [0; 1], ( ) 0, [0; 1]. x f x x ∈ = ∉ (7.1) Функция распределения вероятностей имеет вид 0, 0, ( ) ( ') ' , 0 1, 1, 1. x x F x f x dx x x x −∞ < = = ≤ ≤ > ∫ Числовые характеристики. • математическое ожидание 0,5 m = ; • дисперсия ( ) 1/12 D x = ; • коэффициент асимметрии 0 = A ; коэффициент эксцесса 6 / 5 E = − ; Будем называть x − случайным числом, а ε − случайной цифрой. Установим связь между x и . ε Представим число x в виде бесконечной десятичной дроби 94 1 10 k i k x ε ∞ − = = ⋅ ∑ (7.2) Справедлива следующая теорема: десятичные цифры 1 2 , ,..., k ε ε ε случайного числа x представляют собой независимые случайные цифры. Наоборот, если 1 2 , ,..., k ε ε ε независимые случайные цифры, то формула (7.2) определяет случайное число. Замечание. В вычислениях всегда используют числа с конечным числом десятичных знаков, поэтому случайные числа α заменяют на случайные конечные дроби 1 10 n k i k x ε − = = ⋅ ∑ 7.2. Псевдослучайные числа Пригодность случайных чисел определяется не процессом их получения, а тем, что они должны обладать интересующими нас свойствами независимых, равномерно распределенных СВ. Определение. Последовательность чисел { } i x , которые вычисляются по какой-либо заданной формуле и могут быть использованы вместо случайных чисел при решении задач численным методом, называются псевдослучайными числами Из сказанного следует, что оказываются тождественными те свойства случайных и псевдослучайных чисел, которые требуются для решения широкого круга задач. По отношению к этим задачам разница между физически генерируемыми случайными числами и псевдослучайными практически отсутствует. К преимуществам псевдослучайных чисел можно отнести: – небольшие затраты машинного времени для их получения; – возможность многократного повторного воспроизведения одной и той же последовательности чисел при необходимости; – большой период повторения; – необходимость однократного тестирования алгоритмов вычисления псевдослучайных чисел. Из последнего утверждения следует, что разрабатываемые датчики случайных чисел необходимо подвергать проверке с помощью специальных тестов, которые должны подтверждать их независимость и равномерность распределения. Важной характеристикой последовательности случайных чисел является ее периодичность. Это означает, что имеется некоторый достаточно большой номер L , начиная с которого случайные числа начинают повторяться. Очевидно, что использование при моделировании «большего» отрезка последовательности { } i x , чем период повторения, приведет к бессмысленному повторению испытаний в одних и тех же условиях. 95 7.3. Алгоритмы генераторов псевдослучайных чисел Во всех языках программирования (Pascal, C/C++, Java и т. д.) и в приложениях Excel, MathCad, MathLab и др. есть стандартная функция, возвращающая случайное число. При этом существует возможность повторения одной и той же последовательности случайных чисел, например, в C++, Java. Наиболее распространенные алгоритмы, используемые в генераторах псевдослучайных чисел: 1. Линейный конгруэнтный метод (ЛКМ) − языки Borland С, Visual C++, Java, C++Builder; 2. Алгоритм Вичманна–Хилла (Wichmann–Hill) или AS 183 − языки Prolog, Python (версии 2.2 и предыдущие), Excel; 3. Алгоритм «Виток Мерсенна» (Mersenne Twister) или MT19937−Python (версии 2.3 и последующие); 4. Алгоритм Парка–Миллера; 5. Метод Фибоначчи с запаздыванием ( Subtract-with-borrow Generators SWBG ) −Mathematica, MatLab. Рассмотрим эти алгоритмы. Линейный конгруэнтный метод (ЛКМ) В стандарте ANSI-C имеется функция rand() , выдающая равномерно распределенные числа в интервале от 0 до MAX RAND _ , и связанная с ней функция srand() , выполняющая начальную установку счетчика. Почти все подобные генераторы используют рекуррентную последовательность m c I a I n n mod ) ( 1 + ⋅ = + Здесь 1 + n I равно остатку от деления ( ) c aI n + на m (или другими словами 1 + n I − это наименьший положительный вычет ( ) c aI n + по модулю m ). Число a называется мультипликатором, число c – инкрементом, а число m – модулем. Алгоритм Вичманна–Хилла (Wichmann–Hill или AS183) Псевдослучайные числа вычисляются по формуле 0 , 1 mod 30323 30307 30269 + + = i i i i Z Y X U , где функция 0 , 1 mod возвращает десятичную часть получившейся суммы. Рекуррентные формулы для вычисления i X , i Y и i Z имеют вид: , 30323 mod ) 170 ( , 30307 mod ) 172 ( , 30269 mod ) 171 ( 1 1 1 i i i i i i Z Z Y Y X X × = × = × = + + + где функция mod() возвращает целое число, равное остатку от деления. По сути, этот алгоритм есть линейная комбинация трех конгруэнтных генераторов. При 96 этом требуется задание трёх начальных значений. Алгоритм обладает периодом 12 10 95 , 6 × ( 43 2 ≈ ), что недостаточно для современных нужд. Алгоритм «Виток Мерсенна» (Mersenne Twister или MT19937) Алгоритм разработан в 1997 году японскими учеными Макото Мацумото и Такуджи Нишимура. Обладает огромным периодом 2 19937 −1 (создатели алгоритма доказали это свойство), имеет хорошее быстродействие и проходит все статистические тесты. В приложении 1 приведена реализация алгоритма на языке С. Алгоритм Парка−Миллера (Park, Miller) Самая простая рекуррентная последовательность, которую можно предложить для реализации генератора равномерного распределения, имеет вид: ) mod ) ( ( ) 1 ( m j I a j I × = + Значения констант , 16807 7 5 = = a 31 2 1 2147483647 m = − = были предложены Park, Miller и протестированы в исследованиях Lewis, Goodman, Miller. Прямое использование этого метода возможно на языке ассемблер, но языки высокого уровня могут при этом зафиксировать переполнение. Для обхода этого Scharge предложил метод частичной факторизации модуля. Для этого модуль записывается в виде: r q a + × = mod Если q r < и 1 0 − < < m z , то при этом величины ) mod ( q z a × и ] / [ q z r × всегда лежат в интервале 0,..., m-1. Для вычисления ) mod ( q z a × используется алгоритм: t = a(z mod q)-r[z/q] если t<0, то t += m. (a*z)(mod m)=t. В случае констант Парка−Миллера можно использовать значения 12773 = q и 2836 = r Если требуется число вызовов, превышающее по порядку 10 8 , то для этого случая L'Ecuyer рекомендует комбинировать две последовательностей с близкими, но отличающимися константами. В его исследованиях хороший результат был получен для значений: m1 = 2147483563, a1 = 40014, q1 = 53668, r1 = 12211; m2 = 2147483399, a2 = 40692, q2 = 52774, r2 = 3791. При этом для современных компьютеров период повторения генерируемой последовательности оценивается по порядку примерно как 10 18 97 Метод Фибоначчи с запаздыванием Статистические свойства чисел, генерируемых линейным конгруэнтным алгоритмом, делают невозможным их использование для решения задач, чувствительных к качеству случайных чисел. В связи с этим линейный конгруэнтный алгоритм постепенно потерял свою популярность. Его место заняло семейство фибоначчиевых алгоритмов, позволяющих получать более качественные последовательности псевдослучайных чисел. В англоязычной литературе фибоначчиевы датчики называют обычно «Subtract-with-borrow Generators» (SWBG). Один из фибоначчиевых датчиков основан на следующей итеративной формуле: − < − + − − − − ≥ − − − − = ), ( ) ( если , 1 ) ( ) ( ), ( ) ( если ), ( ) ( ) ( b k X a k X b k X a k X b k X a k X b k X a k X k X где ) (k X − вещественные числа из интервала [0; 1], a, b − целые положительные числа, называемые лагами. Для работы фибоначчиеву датчику требуется знать максимальные значения предыдущих сгенерированных случайных чисел. При программной реализации для хранения полученных чисел используется конечная циклическая очередь на базе массива. Для старта фибоначчиевому датчику требуется ) , max( b a случайных чисел, которые могут быть сгенерированы простым конгруэнтным датчиком. Лаги a и b не следует выбирать произвольно. Рекомендуются следующие пары значений лагов: a = 55, b = 24; a = 17, b = 5; a = 97, b = 33. Качество получаемых случайных чисел зависит от значения константы a, чем оно больше, тем выше размерность пространства, в котором сохраняется равномерность случайных векторов, образованных из полученных случайных чисел. В то же время с увеличением величины константы a возрастает объём используемой алгоритмом памяти. Получаемые случайные числа обладают хорошими статистическими свойствами, причём все биты случайного числа равнозначны по статистическим свойствам. Период фибоначчиева датчика может быть оценен по следующей формуле: ( ) max( , ) 2 1 2 a b j T = − × , где j – число битов в мантиссе вещественного числа. Функция random() в различных приложениях 1. Pascal − random(x) возвращает псевдослучайное число типа word из интервала [0; x]. Переменная x должна быть типа word; 2. Borland C, Microsoft Visual C++ 6.0, Microsoft Visual Studio 2005 − rand(void) возвращает псевдослучайное целое типа int число из интервала [0; 32767]; 3. MathCad − rnd(x) возвращает псевдослучайное число из интервала [0; x]; 4. Excel − (СЛЧИС() – для русской версии) возвращает псевдослучайное число из интервала [0; 1,0]; 5. Mathematica − Random[] возвращает псевдослучайное число типа Real из интервала [0; 1,0]; 98 6. MatLab − rand возвращает случайное псевдослучайное число из интервала [0; 1,0]; 7. Python − import whrandom; print whrandom.random( ) – возвращает псевдослучайное число из интервала [0; 1,0]; 8. Java − Math.random( ) возвращает псевдослучайное число типа double из интервала [0; 1,0]. 7.4 Оценка закона распределения последовательности псевдослучайных чисел Для оценки закона распределения последовательности псевдослучайных чисел, воспользуемся модифицированным критерием Колмогорова- Смирнова. Для этого полученную последовательность случайных чисел расположим в порядке возрастания 1 2 N x x x ≤ ≤ ≤ Рассчитаем следующую статистику 0, 4 0,68 0, 2 D D N N N = + ⋅ + + , (7.3) где ( ) max , D D D + − = ; max ; max 1 1 i i i i i i D x D x n n + − = − = − + + (7.4) Если вычисленная значение статистики не превосходят критического, т.е., если выполняется неравенство D D α ≤ , то равномерность распределения случайных чисел подтверждается. Критическое значение статистики приведено в табл. 7.1. Таблица 7.1 Критическое значение D α критерия равномерности. Уровень значимости α 0,01 0,025 0,05 0,1 0,15 D α 1,628 1,480 1,358 1,224 1,138 7.5. Лабораторное задание по работе №7 1. Для датчика псевдослучайных чисел в приложениях MathCad вычислить оценки среднего, дисперсии, коэффициента асимметрии и коэффициентом эксцесса по формулам 2 1 1 1 1 , ( ) 1 N N i i i i m x D x m N N = = = = − − ∑ ∑ ; 3 3 1 / 2 1 ( ) 1 ( ) N i i A D x m N = = − − ∑ , 4 2 1 1 ( ) 3 1 N i i E D x m N = − − = − ∑ 99 Сравнить полученные оценки с точными значениями математического ожидания 0,5 m = , дисперсии ( ) 1/12 D x = , коэффициента асимметрии 0 = A и коэффициентом эксцесса 6 / 5 E = − Построить гистограмму эмпирической плотности распределения и функции распределения. Объём выборки взять из таблицы 2. 1.1. Проверить гипотезу о равномерности распределения при 0.05 α = 1.2. Получить две выборки объёма N и рассчитать корреляцию. Указания.Для построения гистограммы по выборке объемом n необходимо весь интервал значений случайной величины min max [ , ] x x разбить на m подынтервалов точками i z [ , ), 0,1, 2, ..., i i z z h i m + = , 1 b a h m m − = = (7.13) и вычислить количество i n случайных величин i x , попавших в i -й интервал. Гистограммой относительных частот называется система прямоугольников, каждый из которых имеет основание шириной 1 i i i h z z h + = − = и высоту , 1, 2, ..., i i y i m h ω = = (7.14) где i i n n ω = – площадь i -го прямоугольника . Таблица 7.2. Объём выборки для тестирования датчиков № варианта 1 2 3 4 5 6 7 8 9 10 Объём выборки N 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 2. Для последовательности чисел ( ) 1 0 1 arccos cos(10 ) ; 0,1,..., ; 10; 1. n i i x x i N n x π + = ⋅ = = = Рассчитать эмпирическое среднее, дисперсию, коэффициент асимметрии и коэффициент эксцесса и сравнить с точным значением среднего, дисперсии, коэффициента асимметрии и коэффициентом эксцесса случайной величины с равномерным распределением. Построить гистограмму. Объём выборки взять из таблицы 7. 2.1. Проверить гипотезу о равномерности распределения при 0.05 α = 2.2. Получить две выборки объёма N и рассчитать корреляцию. 100 3. Для последовательности чисел ( ) 1 0 1 arccos cos(10 ) ; 0,1,..., ; 10; 1. n i i x x i N n x π + = ⋅ = = = Рассчитать эмпирическое среднее, дисперсию, коэффициент асимметрии и коэффициент эксцесса и сравнить с точным значением среднего, дисперсии, коэффициента асимметрии и коэффициентом эксцесса случайной величины с равномерным распределением. Построить гистограмму. Объём выборки взять из таблицы 7. 3.1. Проверить гипотезу о равномерности распределения при 0.05 α = 3.1. Получить две выборки объёма N и рассчитать корреляцию. 4. Написать программу линейного конгруэнтного датчика в пакете Mathcad и провести исследование датчика |