Статистическое моделирование
Статистическое определение выходных результатов как основная цель статистического моделирования. Табличные и алгоритмические генераторы случайных чисел. Моделирование случайного события. Моделирование случайной величины с заданным законом распределения.
Рубрика | Математика |
Вид | курс лекций |
Язык | русский |
Дата добавления | 16.04.2013 |
Размер файла | 862,2 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
Размещено на http://www.allbest.ru/
Содержание
- Лекция 1. Статистическое моделирование
- Метод Монте-Карло
- Лекция 2. Генераторы случайных чисел
- Табличные ГСЧ
- Алгоритмические ГСЧ
- Метод серединных квадратов
- Метод серединных произведений
- Метод перемешивания
- Линейный конгруэнтный метод
- Проверка качества работы генератора
- Проверки на равномерность распределения
- Проверки на статистическую независимость
- Лекция 3. Моделирование случайного события
- Моделирование полной группы несовместных событий
- Лекция 4. Моделирование случайной величины с заданным законом распределения
- Метод ступенчатой аппроксимации
- Метод усечения
- Метод взятия обратной функции
Лекция 1. Статистическое моделирование
Статистическое моделирование - базовый метод моделирования, заключающийся в том, что модель испытывается множеством случайных сигналов с заданной плотностью вероятности. Целью является статистическое определение выходных результатов. В основе статистического моделирования лежит метод Монте-Карло. Напомним, что имитацию используют тогда, когда другие методы применить невозможно.
Метод Монте-Карло
Рассмотрим метод Монте-Карло на примере вычисления интеграла, значение которого аналитическим способом найти не удается.
Задача 1. Найти значение интеграла:
На рис. 1.1 представлен график функции f (x). Вычислить значение интеграла этой функции - значит, найти площадь под этим графиком.
Рис. 1.1 Определение значения интеграла методом Монте-Карло
Ограничиваем кривую сверху, справа и слева. Случайным образом распределяем точки в прямоугольнике поиска. Обозначим через N1 количество точек, принятых для испытаний (то есть попавших в прямоугольник, эти точки изображены на рис. 1.1 красным и синим цветом), и через N2 - количество точек под кривой, то есть попавших в закрашенную площадь под функцией (эти точки изображены на рис. 1.1 красным цветом). Тогда естественно предположить, что количество точек, попавших под кривую по отношению к общему числу точек пропорционально площади под кривой (величине интеграла) по отношению к площади испытуемого прямоугольника. Математически это можно выразить так:
Рассуждения эти, конечно, статистические и тем более верны, чем большее число испытуемых точек мы возьмем.
Фрагмент алгоритма метода Монте-Карло в виде блок-схемы выглядит так, как показано на рис. 1.2
Рис. 1.2 Фрагмент алгоритма реализации метода Монте-Карло
Значения r1 и r2 на рис. 1.2 являются равномерно распределенными случайными числами из интервалов (x1; x2) и (c1; c2) соответственно.
Метод Монте-Карло чрезвычайно эффективен, прост, но необходим "хороший" генератор случайных чисел. Вторая проблема применения метода заключается в определении объема выборки, то есть количества точек, необходимых для обеспечения решения с заданной точностью. Эксперименты показывают: чтобы увеличить точность в 10 раз, объем выборки нужно увеличить в 100 раз; то есть точность примерно пропорциональна корню квадратному из объема выборки:
Схема использования метода Монте-Карло при исследовании систем со случайными параметрами
Построив модель системы со случайными параметрами, на ее вход подают входные сигналы от генератора случайных чисел (ГСЧ), как показано на рис. 1.3 ГСЧ устроен так, что он выдает равномерно распределенные случайные числа rрр из интервала [0; 1]. Так как одни события могут быть более вероятными, другие - менее вероятными, то равномерно распределенные случайные числа от генератора подают на преобразователь закона случайных чисел (ПЗСЧ), который преобразует их в заданный пользователем закон распределения вероятности, например, в нормальный или экспоненциальный закон. Эти преобразованные случайные числа x подают на вход модели. Модель отрабатывает входной сигнал x по некоторому закону y = ц (x) и получает выходной сигнал y, который также является случайным.
статистическое моделирование случайная величина
Рис. 1.3 Общая схема метода статистического моделирования
В блоке накопления статистики (БНСтат) установлены фильтры и счетчики. Фильтр (некоторое логическое условие) определяет по значению y, реализовалось ли в конкретном опыте некоторое событие (выполнилось условие, f = 1) или нет (условие не выполнилось, f = 0). Если событие реализовалось, то счетчик события увеличивается на единицу. Если событие не реализовалось, то значение счетчика не меняется. Если требуется следить за несколькими разными типами событий, то для статистического моделирования понадобится несколько фильтров и счетчиков Ni. Всегда ведется счетчик количества экспериментов - N.
Далее отношение Ni к N, рассчитываемое в блоке вычисления статистических характеристик (БВСХ) по методу Монте-Карло, дает оценку вероятности pi появления события i, то есть указывает на частоту его выпадения в серии из N опытов. Это позволяет сделать выводы о статистических свойствах моделируемого объекта.
Например, событие A совершилось в результате проведенных 200 экспериментов 50 раз. Это означает, согласно методу Монте-Карло, что вероятность совершения события равна: pA = 50/200 = 0.25. Вероятность того, что событие не совершится, равна, соответственно, 1 - 0.25 = 0.75.
Обратите внимание: когда говорят о вероятности, полученной экспериментально, то ее называют частостью; слово вероятность употребляют, когда хотят подчеркнуть, что речь идет о теоретическом понятии.
При большом количестве опытов N частота появления события, полученная экспериментальным путем, стремится к значению теоретической вероятности появления события.
В блоке оценки достоверности (БОД) анализируют степень достоверности статистических экспериментальных данных, снятых с модели (принимая во внимание точность результата е, заданную пользователем) и определяют необходимое для этого количество статистических испытаний. Если колебания значений частоты появления событий относительно теоретической вероятности меньше заданной точности, то экспериментальную частоту принимают в качестве ответа, иначе генерацию случайных входных воздействий продолжают, и процесс моделирования повторяется. При малом числе испытаний результат может оказаться недостоверным. Но чем более испытаний, тем точнее ответ, согласно центральной предельной теореме.
Заметим, что оценивание ведут по худшей из частот. Это обеспечивает достоверный результат сразу по всем снимаемым характеристикам модели.
Пример 1. Решим простую задачу. Какова вероятность выпадения монеты орлом кверху при падении ее с высоты случайным образом?
Начнем подбрасывать монетку и фиксировать результаты каждого броска (см. табл. 1.1).
Таблица 1.1.
Результаты испытаний бросания монеты
Будем подсчитывать частость выпадения орла как отношение количества случаев выпадения орла к общему числу наблюдений. Посмотрите в табл. 1.1 случаи для N = 1, N = 2, N = 3 - сначала значения частости нельзя назвать достоверными. Попробуем построить график зависимости Pо от N - и посмотрим, как меняется частость выпадения орла в зависимости от количества проведенных опытов. Разумеется, при различных экспериментах будут получаться разные таблицы и, следовательно, разные графики. На рис. 1.4 показан один из вариантов.
Рис. 1.4 Экспериментальная зависимость частости появления случайного события от количества наблюдений и ее стремление к теоретической вероятности
Сделаем некоторые выводы.
1. Видно, что при малых значениях N, например, N = 1, N = 2, N = 3 ответу вообще доверять нельзя. Например, Pо = 0 при N = 1, то есть вероятность выпадения орла при одном броске равна нулю! Хотя всем хорошо известно, что это не так. То есть пока мы получили очень грубый ответ. Однако, посмотрите на график: в процессе накопления информации ответ медленно, но верно приближается к правильному (он выделен пунктирной линией). К счастью, в данном конкретном случае правильный ответ нам известен: в идеале, вероятность выпадения орла равна 0.5 (в других, более сложных задачах, ответ нам, конечно, будет неизвестен). Допустим, что ответ нам надо знать с точностью е = 0.1. Проведем две параллельные линии, отстоящие от правильного ответа 0.5 на расстояние 0.1 (см. рис. 1.4). Ширина образовавшегося коридора будет равна 0.2. Как только кривая Pо (N) войдет в этот коридор так, что уже никогда его не покинет, можно остановиться и посмотреть, для какого значения N это произошло. Это и есть экспериментально вычисленное критическое значение необходимого количества опытов Nкрэ для определения ответа с точностью е = 0.1; е-окрестность в наших рассуждениях играет роль своеобразной трубки точности. Заметьте, что ответы Pо (91), Pо (92) и так далее уже не меняют сильно своих значений (см. рис. 1.4); по крайней мере, у них не изменяется первая цифра после запятой, которой мы обязаны доверять по условиям задачи.
2. Причиной такого поведения кривой является действие центральной предельной теоремы. Пока здесь мы сформулируем ее в самом простом варианте "Сумма случайных величин есть величина неслучайная". Мы использовали среднюю величину Pо, которая несет в себе информацию о сумме опытов, и поэтому постепенно эта величина становится все более достоверной.
3. Если проделать еще раз этот опыт сначала, то, конечно, его результатом будет другой вид случайной кривой. И ответ будет другим, хотя примерно таким же. Проведем целую серию таких экспериментов (см. рис. 1.5). Такая серия называется ансамблем реализаций. Какому же ответу в итоге следует верить? Ведь они, хоть и являются близкими, все же разнятся. На практике поступают по-разному. Первый вариант - вычислить среднее значение ответов за несколько реализаций (см. табл. 1.2).
Рис. 1.5 Экспериментально снятый ансамбль случайных зависимостей частости появления случайного события от количества наблюдений
Мы поставили несколько экспериментов и определяли каждый раз, сколько необходимо было сделать опытов, то есть Nкрэ. Было проделано 10 экспериментов, результаты которых были сведены в табл. 1.2 По результатам 10-ти экспериментов было вычислено среднее значение Nкрэ.
Таблица 1.2.
Экспериментальные данные необходимого количества бросков монеты для достижения точности е = 0.1 при вычислении вероятности выпадения орла
Таким образом, проведя 10 реализаций разной длины, мы определили, что достаточно в среднем было сделать 1 реализацию длиной в 94 броска монеты.
Еще один важный факт. Внимательно рассмотрите график на рис.21.5 На нем нарисовано 100 реализаций - 100 красных линий. Отметьте на нем абсциссу N = 94 вертикальной чертой. Есть какой-то процент красных линий, которые не успели пересечь е-окрестность, то есть (Pэксп - е ? Pтеор ? Pэксп + е), и войти в коридор точности до момента N = 94. Обратите внимание, таких линий 5. Это значит, что 95 из 100, то есть 95%, линий достоверно вошли в обозначенный интервал.
Таким образом, проведя 100 реализаций, мы добились примерно 95% -ного доверия к полученной экспериментально величине вероятности выпадения орла, определив ее с точностью 0.1.
Для сравнения полученного результата вычислим теоретическое значение Nкрт теоретически. Однако для этого придется ввести понятие доверительной вероятности QF, которая показывает, насколько мы готовы верить ответу.
Например, при QF = 0.95 мы готовы верить ответу в 95% случаев из 100. Имеет вид: Nкрт = k (QF) · p · (1 - p) /е2, где k (QF) - коэффициент Лапласа, p - вероятность выпадения орла, е - точность (доверительный интервал). В табл. 1.3 показаны значения теоретической величины количества необходимых опытов при разных QF (для точности е = 0.1 и вероятности p = 0.5).
Таблица 1.3.
Теоретический расчет необходимого количества бросков монеты для достижения точности е = 0.1 при вычислении вероятности выпадения орла
Как видите, полученная нами оценка длины реализации, равная 94 опытам очень близка к теоретической, равной 96. Некоторое несовпадение объясняется тем, что, видимо, 10 реализаций недостаточно для точного вычисления Nкрэ. Если вы решите, что вам нужен результат, которому следует доверять больше, то измените значение доверительной вероятности. Например, теория говорит нам, что если опытов будет 167, то всего 1-2 линии из ансамбля не войдут в предложенную трубку точности. Но имейте в виду, количество экспериментов с ростом точности и достоверности растет очень быстро.
Второй вариант, используемый на практике - провести одну реализацию и увеличить полученное для нее Nкрэ в 2 раза. Это считают хорошей гарантией точности ответа (см. рис. 1.6).
Рис. 1.6. Иллюстрация экспериментального определения Nкрэ по правилу "умножь на два"
Если присмотреться к ансамблю случайных реализаций, то можно обнаружить, что сходимость частости к значению теоретической вероятности происходит по кривой, соответствующей обратной квадратичной зависимости от числа экспериментов (см. рис. 1.7).
Рис. 1.7 Иллюстрация скорости схождения экспериментально получаемой частости к теоретической вероятности
Это действительно так получается и теоретически. Если изменять задаваемую точность е и исследовать количество экспериментов, требуемых для обеспечения каждой из них, то получится табл. 1.4
Таблица 1.4.
Теоретическая зависимость количества экспериментов, необходимых для обеспечения заданной точности при QF = 0.95
Построим по табл. 1.4 график зависимости Nкрт (е) (см. рис. 1.8).
Рис. 1.8 Зависимость числа экспериментов, требуемых для достижения заданной точности е при фиксированном QF = 0.95
Итак, рассмотренные графики подтверждают приведенную выше оценку:
Заметим, что оценок точности может быть несколько.
Пример 2. Нахождение площади фигуры методом Монте-Карло. Определите методом Монте-Карло площадь пятиугольника с координатами углов (0, 0), (0,10), (5, 20), (10,10), (7, 0).
Нарисуем в двухмерных координатах заданный пятиугольник, вписав его в прямоугольник, чья площадь, как нетрудно догадаться, составляет (10 - 0) · (20 - 0) = 200 (см. рис. 1.9).
Рис. 1.9 Иллюстрация к решению задачи о площади фигуры методом Монте-Карло
Используем таблицу случайных чисел для генерации пар чисел R, G, равномерно распределенных в интервале от 0 до 1. Число R будет имитировать координату X (0 ? X ? 10), следовательно, X = 10 · R. Число G будет имитировать координату Y (0 ? Y ? 20), следовательно, Y = 20 · G. Сгенерируем по 10 чисел R и G и отобразим 10 точек (X; Y) на рис. 1.9 и в табл. 1.5
Таблица 1.5.
Решение задачи методом Монте-Карло
Статистическая гипотеза заключается в том, что количество точек, попавших в контур фигуры, пропорционально площади фигуры: 6: 10 = S: 200. То есть, по формуле метода Монте-Карло, получаем, что площадь S пятиугольника равна: 200 · 6/10 = 120.
Проследим, как менялась величина S от опыта к опыту (см. табл. 1.6).
Таблица 1.6.
Оценка точности ответа
Количество испытаний N |
Оценка вероятности попадания случайной точки в испытуемую область |
Оценка площади S методом Монте-Карло |
|
1 |
1/1 = 1.00 |
200 |
|
2 |
1/2 = 0.50 |
100 |
|
3 |
2/3 = 0.67 |
133 |
|
4 |
3/4 = 0.75 |
150 |
|
5 |
3/5 = 0.60 |
120 |
|
6 |
4/6 = 0.67 |
133 |
|
7 |
5/7 = 0.71 |
143 |
|
8 |
5/8 = 0.63 |
125 |
|
9 |
6/9 = 0.67 |
133 |
|
10 |
6/10 = 0.60 |
120 |
Поскольку в ответе все еще меняется значение второго разряда, то возможная неточность составляет пока больше 10%. Точность расчета может быть увеличена с ростом числа испытаний (см. рис. 1.10).
Рис. 1.10 Иллюстрация процесса сходимости определяемого экспериментально ответа к теоретическому результату
Лекция 2. Генераторы случайных чисел
В основе метода Монте-Карло (см. Лекцию 1. Статистическое моделирование) лежит генерация случайных чисел, которые должны быть равномерно распределены в интервале (0;1).
Если генератор выдает числа, смещенные в какую-то часть интервала (одни числа выпадают чаще других), то результат решения задачи, решаемой статистическим методом, может оказаться неверным. Поэтому проблема использования хорошего генератора действительно случайных и действительно равномерно распределенных чисел стоит очень остро.
Математическое ожидание mr и дисперсия Dr такой последовательности, состоящей из n случайных чисел ri, должны быть следующими (если это действительно равномерно распределенные случайные числа в интервале от 0 до 1):
Если пользователю потребуется, чтобы случайное число x находилось в интервале (a; b), отличном от (0;
1), нужно воспользоваться формулой x = a + (b - a) · r, где r - случайное число из интервала (0;
1). Законность данного преобразования демонстрируется на рис. 2.1
Рис. 2.1 Схема перевода числа из интервала (0;
1) в интервал (a; b)
Теперь x - случайное число, равномерно распределенное в диапазоне от a до b.
За эталон генератора случайных чисел (ГСЧ) принят такой генератор, который порождает последовательность случайных чисел с равномерным законом распределения в интервале (0;
1). За одно обращение данный генератор возвращает одно случайное число. Если наблюдать такой ГСЧ достаточно длительное время, то окажется, что, например, в каждый из десяти интервалов (0; 0.1), (0.1; 0.2), (0.2; 0.3), …, (0.9;
1) попадет практически одинаковое количество случайных чисел - то есть они будут распределены равномерно по всему интервалу (0;
1). Если изобразить на графике k = 10 интервалов и частоты Ni попаданий в них, то получится экспериментальная кривая плотности распределения случайных чисел (см. рис. 2.2).
Рис. 2.2 Частотная диаграмма выпадения случайных чисел, порождаемых реальным генератором
Заметим, что в идеале кривая плотности распределения случайных чисел выглядела бы так, как показано на рис. 2.3. То есть в идеальном случае в каждый интервал попадает одинаковое число точек: Ni = N/k, где N - общее число точек, k - количество интервалов, i = 1, …, k.
Рис. 2.3 Частотная диаграмма выпадения случайных чисел, порождаемых идеальным генератором теоретически
Следует помнить, что генерация произвольного случайного числа состоит из двух этапов:
· генерация нормализованного случайного числа (то есть равномерно распределенного от 0 до 1);
· преобразование нормализованных случайных чисел ri в случайные числа xi, которые распределены по необходимому пользователю (произвольному) закону распределения или в необходимом интервале.
Генераторы случайных чисел по способу получения чисел делятся на:
· физические;
· табличные;
· алгоритмические.
Физические ГСЧ
Примером физических ГСЧ могут служить: монета ("орел" - 1, "решка" - 0); игральные кости; поделенный на секторы с цифрами барабан со стрелкой; аппаратурный генератор шума (ГШ), в качестве которого используют шумящее тепловое устройство, например, транзистор (рис. 2.4-2.5).
Рис. 2.4 Схема аппаратного метода генерации случайных чисел
Рис. 2.5 Диаграмма получения случайных чисел аппаратным методом
Задача "Генерация случайных чисел при помощи монеты"
Сгенерируйте случайное трехразрядное число, распределенное по равномерному закону в интервале от 0 до 1, с помощью монеты. Точность - три знака после запятой.
Первый способ решения задачи Подбросьте монету 9 раз, и если монета упала решкой, то запишите "0", если орлом, то "1". Итак, допустим, что в результате эксперимента получили случайную последовательность 100110100.
Начертите интервал от 0 до 1. Считывая числа в последовательности слева направо, разбивайте интервал пополам и выбирайте каждый раз одну из частей очередного интервала (если выпал 0, то левую, если выпала 1, то правую). Таким образом, можно добраться до любой точки интервала, сколь угодно точно.
Итак, 1: интервал [0; 1] делится пополам - [0; 0.5] и [0.5; 1], - выбирается правая половина, интервал сужается: [0.5; 1]. Следующее число, 0: интервал [0.5; 1] делится пополам - [0.5; 0.75] и [0.75; 1], - выбирается левая половина [0.5; 0.75], интервал сужается: [0.5; 0.75]. Следующее число, 0: интервал [0.5; 0.75] делится пополам - [0.5; 0.625] и [0.625; 0.75], - выбирается левая половина [0.5; 0.625], интервал сужается: [0.5; 0.625]. Следующее число, 1: интервал [0.5; 0.625] делится пополам - [0.5; 0.5625] и [0.5625; 0.625], - выбирается правая половина [0.5625; 0.6250], интервал сужается: [0.5625; 0.6250].
По условию точности задачи решение найдено: им является любое число из интервала [0.5625; 0.6250], например, 0.625.
В принципе, если подходить строго, то деление интервалов нужно продолжить до тех пор, пока левая и правая границы найденного интервала не СОВПАДУТ между собой с точностью до третьего знака после запятой. То есть с позиций точности сгенерированное число уже не будет отличимо от любого числа из интервала, в котором оно находится.
Второй способ решения задачи Разобьем полученную двоичную последовательность 100110100 на триады: 100, 110, 100. После перевода этих двоичных чисел в десятичные получаем: 4, 6,4. Подставив спереди "0.", получим: 0.464. Таким методом могут получаться только числа от 0.000 до 0.777 (так как максимум, что можно "выжать" из трех двоичных разрядов - это 1112 = 78) - то есть, по сути, эти числа представлены в восьмеричной системе счисления. Для перевода восьмеричного числа в десятичное представление выполним:
0.4648 = 4 · 8-1 + 6 · 8-2 + 4 · 8-3 = 0.601562510 = 0.60210.
Итак, искомое число равно: 0.602.
Табличные ГСЧ
Табличные ГСЧ в качестве источника случайных чисел используют специальным образом составленные таблицы, содержащие проверенные некоррелированные, то есть никак не зависящие друг от друга, цифры. В табл.22.1 приведен небольшой фрагмент такой таблицы. Обходя таблицу слева направо сверху вниз, можно получать равномерно распределенные от 0 до 1 случайные числа с нужным числом знаков после запятой (в нашем примере мы используем для каждого числа по три знака). Так как цифры в таблице не зависят друг от друга, то таблицу можно обходить разными способами, например, сверху вниз, или справа налево, или, скажем, можно выбирать цифры, находящиеся на четных позициях.
Таблица 2.1.
Случайные цифры. Равномерно распределенные от 0 до 1 случайные числа
Достоинство данного метода в том, что он дает действительно случайные числа, так как таблица содержит проверенные некоррелированные цифры. Недостатки метода: для хранения большого количества цифр требуется много памяти; большие трудности порождения и проверки такого рода таблиц, повторы при использовании таблицы уже не гарантируют случайности числовой последовательности, а значит, и надежности результата.
Здесь находится таблица, содержащая 500 абсолютно случайных проверенных чисел (взято из книги И.Г. Венецкого, В.И. Венецкой "Основные математико-статистические понятия и формулы в экономическом анализе").
Алгоритмические ГСЧ
Числа, генерируемые с помощью этих ГСЧ, всегда являются псевдослучайными (или квазислучайными), то есть каждое последующее сгенерированное число зависит от предыдущего:
ri + 1 = f (ri).
Последовательности, составленные из таких чисел, образуют петли, то есть обязательно существует цикл, повторяющийся бесконечное число раз. Повторяющиеся циклы называются периодами.
Достоинством данных ГСЧ является быстродействие; генераторы практически не требуют ресурсов памяти, компактны. Недостатки: числа нельзя в полной мере назвать случайными, поскольку между ними имеется зависимость, а также наличие периодов в последовательности квазислучайных чисел.
Рассмотрим несколько алгоритмических методов получения ГСЧ:
· метод серединных квадратов;
· метод серединных произведений;
· метод перемешивания;
· линейный конгруэнтный метод.
Метод серединных квадратов
Имеется некоторое четырехзначное число R0. Это число возводится в квадрат и заносится в R1. Далее из R1 берется середина (четыре средних цифры) - новое случайное число - и записывается в R0. Затем процедура повторяется (см. рис. 2.6). Отметим, что на самом деле в качестве случайного числа необходимо брать не ghij, а 0. ghij - с приписанным слева нулем и десятичной точкой. Этот факт отражен как на рис. 2.6, так и на последующих подобных рисунках.
Рис. 2.6 Схема метода серединных квадратов
Недостатки метода:
1) если на некоторой итерации число R0 станет равным нулю, то генератор вырождается, поэтому важен правильный выбор начального значения R0;
2) генератор будет повторять последовательность через Mn шагов (в лучшем случае), где n - разрядность числа R0, M - основание системы счисления.
Для примера на рис. 2.6: если число R0 будет представлено в двоичной системе счисления, то последовательность псевдослучайных чисел повторится через 24 = 16 шагов. Заметим, что повторение последовательности может произойти и раньше, если начальное число будет выбрано неудачно.
Описанный выше способ был предложен Джоном фон Нейманом и относится к 1946 году. Поскольку этот способ оказался ненадежным, от него очень быстро отказались.
Метод серединных произведений
Число R0 умножается на R1, из полученного результата R2 извлекается середина R2* (это очередное случайное число) и умножается на R1. По этой схеме вычисляются все последующие случайные числа (см. рис. 2.7).
Рис. 2.7 Схема метода серединных произведений
Метод перемешивания
В методе перемешивания используются операции циклического сдвига содержимого ячейки влево и вправо. Идея метода состоит в следующем. Пусть в ячейке хранится начальное число R0. Циклически сдвигая содержимое ячейки влево на 1/4 длины ячейки, получаем новое число R0*. Точно так же, циклически сдвигая содержимое ячейки R0 вправо на 1/4 длины ячейки, получаем второе число R0**. Сумма чисел R0* и R0** дает новое случайное число R1. Далее R1 заносится в R0, и вся последовательность операций повторяется (см. рис. 2.8).
Рис. 2.8 Схема метода перемешивания
Обратите внимание, что число, полученное в результате суммирования R0* и R0**, может не уместиться полностью в ячейке R1. В этом случае от полученного числа должны быть отброшены лишние разряды. Поясним это для рис.22.8, где все ячейки представлены восемью двоичными разрядами. Пусть R0* = 100100012 = 14510, R0** = 101000012 = 16110, тогда R0* + R0** = 1001100102 = 30610. Как видим, число 306 занимает 9 разрядов (в двоичной системе счисления), а ячейка R1 (как и R0) может вместить в себя максимум 8 разрядов. Поэтому перед занесением значения в R1 необходимо убрать один "лишний", крайний левый бит из числа 306, в результате чего в R1 пойдет уже не 306, а 001100102 = 5010. Также заметим, что в таких языках, как Паскаль, "урезание" лишних битов при переполнении ячейки производится автоматически в соответствии с заданным типом переменной.
Линейный конгруэнтный метод
Линейный конгруэнтный метод является одной из простейших и наиболее употребительных в настоящее время процедур, имитирующих случайные числа. В этом методе используется операция mod (x, y), возвращающая остаток от деления первого аргумента на второй. Каждое последующее случайное число рассчитывается на основе предыдущего случайного числа по следующей формуле:
ri + 1 = mod (k · ri + b, M).
M - модуль (0 < M);
k - множитель (0 ? k < M);
b - приращение (0 ? b < M);
r0 - начальное значение (0 ? r0 < M).
Последовательность случайных чисел, полученных с помощью данной формулы, называется линейной конгруэнтной последовательностью. Многие авторы называют линейную конгруэнтную последовательность при b = 0 мультипликативным конгруэнтным методом, а при b ? 0 - смешанным конгруэнтным методом.
Для качественного генератора требуется подобрать подходящие коэффициенты. Необходимо, чтобы число M было довольно большим, так как период не может иметь больше M элементов. С другой стороны, деление, использующееся в этом методе, является довольно медленной операцией, поэтому для двоичной вычислительной машины логичным будет выбор M = 2N, поскольку в этом случае нахождение остатка от деления сводится внутри ЭВМ к двоичной логической операции "AND". Также широко распространен выбор наибольшего простого числа M, меньшего, чем 2N: в специальной литературе доказывается, что в этом случае младшие разряды получаемого случайного числа ri + 1 ведут себя так же случайно, как и старшие, что положительно сказывается на всей последовательности случайных чисел в целом. В качестве примера можно привести одно из чисел Мерсенна, равное 231 - 1, и таким образом, M = 231 - 1.
Одним из требований к линейным конгруэнтным последовательностям является как можно большая длина периода. Длина периода зависит от значений M, k и b. Теорема, которую мы приведем ниже, позволяет определить, возможно ли достижение периода максимальной длины для конкретных значений M, k и b.
Теорема. Линейная конгруэнтная последовательность, определенная числами M, k, b и r0, имеет период длиной M тогда и только тогда, когда:
· числа b и M взаимно простые;
· k - 1 кратно p для каждого простого p, являющегося делителем M;
· k - 1 кратно 4, если M кратно 4.
Наконец, в заключение рассмотрим пару примеров использования линейного конгруэнтного метода для генерации случайных чисел.
Пример 1
M = 2N k = 3 + 8 · q (или k = 5 + 8 · q) b = 0 r0 - нечетно
Было установлено, что ряд псевдослучайных чисел, генерируемых на основе данных из примера 1, будет повторяться через каждые M/4 чисел. Число q задается произвольно перед началом вычислений, однако при этом следует иметь в виду, что ряд производит впечатление случайного при больших k (а значит, и q). Результат можно несколько улучшить, если b нечетно и k = 1 + 4 · q - в этом случае ряд будет повторяться через каждые M чисел. После долгих поисков k исследователи остановились на значениях 69069 и 71365.
Пример 2
M = 231 - 1 k = 1 220 703 125 b = 7 r0 = 7
Генератор случайных чисел, использующий данные из примера 2, будет выдавать случайные неповторяющиеся числа с периодом, равным 7 миллионам. Мультипликативный метод генерации псевдослучайных чисел был предложен Д.Г. Лехмером (D. H. Lehmer) в 1949 году.
Проверка качества работы генератора
От качества работы ГСЧ зависит качество работы всей системы и точность результатов. Поэтому случайная последовательность, порождаемая ГСЧ, должна удовлетворять целому ряду критериев.
Осуществляемые проверки бывают двух типов:
· проверки на равномерность распределения;
· проверки на статистическую независимость.
Проверки на равномерность распределения
1) ГСЧ должен выдавать близкие к следующим значения статистических параметров, характерных для равномерного случайного закона:
- математическое ожидание;
- дисперсия;
- среднеквадратичное отклонение.
2) Частотный тест
Частотный тест позволяет выяснить, сколько чисел попало в интервал (mr - уr; mr + уr), то есть (0.5 - 0.2887; 0.5 + 0.2887) или, в конечном итоге, (0.2113; 0.7887). Так как 0.7887 - 0.2113 = 0.5774, заключаем, что в хорошем ГСЧ в этот интервал должно попадать около 57.7% из всех выпавших случайных чисел (см. рис. 2.9).
Рис. 2.9 Частотная диаграмма идеального ГСЧ в случае проверки его на частотный тест
Также необходимо учитывать, что количество чисел, попавших в интервал (0; 0.5), должно быть примерно равно количеству чисел, попавших в интервал (0.5;1).
3) Проверка по критерию "хи-квадрат"
Критерий "хи-квадрат" (ч2-критерий) - это один из самых известных статистических критериев; он является основным методом, используемым в сочетании с другими критериями. Критерий "хи-квадрат" был предложен в 1900 году Карлом Пирсоном. Его замечательная работа рассматривается как фундамент современной математической статистики.
Для нашего случая проверка по критерию "хи-квадрат" позволит узнать, насколько созданный нами реальный ГСЧ близок к эталону ГСЧ, то есть удовлетворяет ли он требованию равномерного распределения или нет.
Частотная диаграмма эталонного ГСЧ представлена на рис. 2.10.
Так как закон распределения эталонного ГСЧ равномерный, то (теоретическая) вероятность pi попадания чисел в i-ый интервал (всего этих интервалов k) равна pi = 1/k.
И, таким образом, в каждый из k интервалов попадет ровно по pi · N чисел (N - общее количество сгенерированных чисел).
Рис. 2.10. Частотная диаграмма эталонного ГСЧ
Реальный ГСЧ будет выдавать числа, распределенные (причем, не обязательно равномерно!) по k интервалам и в каждый интервал попадет по ni чисел (в сумме n1 + n2 + … + nk = N).
Как же нам определить, насколько испытываемый ГСЧ хорош и близок к эталонному?
Вполне логично рассмотреть квадраты разностей между полученным количеством чисел ni и "эталонным" pi · N. Сложим их, и в результате получим:
ч2эксп. = (n1 - p1 · N) 2 + (n2 - p2 · N) 2 + … + (nk - pk · N) 2.
Из этой формулы следует, что чем меньше разность в каждом из слагаемых (а значит, и чем меньше значение ч2эксп.), тем сильнее закон распределения случайных чисел, генерируемых реальным ГСЧ, тяготеет к равномерному.
В предыдущем выражении каждому из слагаемых приписывается одинаковый вес (равный 1), что на самом деле может не соответствовать действительности; поэтому для статистики "хи-квадрат" необходимо провести нормировку каждого i-го слагаемого, поделив его на pi · N:
Наконец, запишем полученное выражение более компактно и упростим его:
Мы получили значение критерия "хи-квадрат" для экспериментальных данных.
В табл. 2.2 приведены теоретические значения "хи-квадрат" (ч2теор.), где н = N - 1 - это число степеней свободы, p - это доверительная вероятность, задаваемая пользователем, который указывает, насколько ГСЧ должен удовлетворять требованиям равномерного распределения, или p - это вероятность того, что экспериментальное значение ч2эксп. будет меньше табулированного (теоретического) ч2теор. или равно ему.
Таблица 22.2.
Некоторые процентные точки ч2-распределения
p = 1% |
p = 5% |
p = 25% |
p = 50% |
p = 75% |
p = 95% |
p = 99% |
||
н = 1 |
0.00016 |
0.00393 |
0.1015 |
0.4549 |
1.323 |
3.841 |
6.635 |
|
н = 2 |
0.02010 |
0.1026 |
0.5754 |
1.386 |
2.773 |
5.991 |
9.210 |
|
н = 3 |
0.1148 |
0.3518 |
1.213 |
2.366 |
4.108 |
7.815 |
11.34 |
|
н = 4 |
0.2971 |
0.7107 |
1.923 |
3.357 |
5.385 |
9.488 |
13.28 |
|
н = 5 |
0.5543 |
1.1455 |
2.675 |
4.351 |
6.626 |
11.07 |
15.09 |
|
н = 6 |
0.8721 |
1.635 |
3.455 |
5.348 |
7.841 |
12.59 |
16.81 |
|
н = 7 |
1.239 |
2.167 |
4.255 |
6.346 |
9.037 |
14.07 |
18.48 |
|
н = 8 |
1.646 |
2.733 |
5.071 |
7.344 |
10.22 |
15.51 |
20.09 |
|
н = 9 |
2.088 |
3.325 |
5.899 |
8.343 |
11.39 |
16.92 |
21.67 |
|
н = 10 |
2.558 |
3.940 |
6.737 |
9.342 |
12.55 |
18.31 |
23.21 |
|
н = 11 |
3.053 |
4.575 |
7.584 |
10.34 |
13.70 |
19.68 |
24.72 |
|
н = 12 |
3.571 |
5.226 |
8.438 |
11.34 |
14.85 |
21.03 |
26.22 |
|
н = 15 |
5.229 |
7.261 |
11.04 |
14.34 |
18.25 |
25.00 |
30.58 |
|
н = 20 |
8.260 |
10.85 |
15.45 |
19.34 |
23.83 |
31.41 |
37.57 |
|
н = 30 |
14.95 |
18.49 |
24.48 |
29.34 |
34.80 |
43.77 |
50.89 |
|
н = 50 |
29.71 |
34.76 |
42.94 |
49.33 |
56.33 |
67.50 |
76.15 |
|
н > 30 |
н + sqrt(2н) · xp + 2/3 · x2p - 2/3 + O(1/sqrt(н)) |
|||||||
xp = |
-2.33 |
-1.64 |
-0.674 |
0.00 |
0.674 |
1.64 |
2.33 |
Приемлемым считают p от 10% до 90%.
Если ч2эксп. много больше ч2теор. (то есть p - велико), то генератор не удовлетворяет требованию равномерного распределения, так как наблюдаемые значения ni слишком далеко уходят от теоретических pi · N и не могут рассматриваться как случайные. Другими словами, устанавливается такой большой доверительный интервал, что ограничения на числа становятся очень нежесткими, требования к числам - слабыми. При этом будет наблюдаться очень большая абсолютная погрешность.
Еще Д. Кнут в своей книге "Искусство программирования" заметил, что иметь ч2эксп. маленьким тоже, в общем-то, нехорошо, хотя это и кажется, на первый взгляд, замечательно с точки зрения равномерности. Действительно, возьмите ряд чисел 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, … - они идеальны с точки зрения равномерности, и ч2эксп. будет практически нулевым, но вряд ли вы их признаете случайными.
Если ч2эксп. много меньше ч2теор. (то есть p - мало), то генератор не удовлетворяет требованию случайного равномерного распределения, так как наблюдаемые значения ni слишком близки к теоретическим pi · N и не могут рассматриваться как случайные.
А вот если ч2эксп. лежит в некотором диапазоне, между двумя значениями ч2теор., которые соответствуют, например, p = 25% и p = 50%, то можно считать, что значения случайных чисел, порождаемые датчиком, вполне являются случайными.
При этом дополнительно надо иметь в виду, что все значения pi · N должны быть достаточно большими, например больше 5 (выяснено эмпирическим путем). Только тогда (при достаточно большой статистической выборке) условия проведения эксперимента можно считать удовлетворительными.
Итак, процедура проверки имеет следующий вид.
1. Диапазон от 0 до 1 разбивается на k равных интервалов.
2. Запускается ГСЧ N раз (N должно быть велико, например, N/k > 5).
3. Определяется количество случайных чисел, попавших в каждый интервал: ni, i = 1, …, k.
4. Вычисляется экспериментальное значение ч2эксп. по следующей формуле:
где pi = 1/k - теоретическая вероятность попадания чисел в k-ый интервал.
5. Путем сравнения экспериментально полученного значения ч2эксп. с теоретическим ч2теор. (из табл. 2.2) делается вывод о пригодности генератора для использования. Для этого: а) входим в табл. 2.2 (строка = количество экспериментов - 1); б) сравниваем вычисленное ч2эксп. с ч2теор., встречающимися в строке. При этом возможно три случая.
Первый случай: ч2эксп. много больше любого ч2теор. в строке - гипотеза о случайности равномерного генератора не выполняется (разброс чисел слишком велик, чтобы быть случайным).
Второй случай: ч2эксп. много меньше любого ч2теор. в строке - гипотеза о случайности равномерного генератора не выполняется (разброс чисел слишком мал, чтобы быть случайным).
Третий случай: ч2эксп. лежит между значениями ч2теор. двух рядом стоящих столбцов - гипотеза о случайности равномерного генератора выполняется с вероятностью p (то есть в p случаях из 100).
Заметим, что чем ближе получается p к значению 50%, тем лучше.
Проверки на статистическую независимость
1) Проверка на частоту появления цифры в последовательности
Рассмотрим пример. Случайное число 0.2463389991 состоит из цифр 2463389991, а число 0.5467766618 состоит из цифр 5467766618. Соединяя последовательности цифр, имеем: 24633899915467766618.
Понятно, что теоретическая вероятность pi выпадения i-ой цифры (от 0 до 9) равна 0.1.
Далее следует вычислить частоту появления каждой цифры в выпавшей экспериментальной последовательности. Например, цифра 1 выпала 2 раза из 20, а цифра 6 выпала 5 раз из 20.
Далее считают оценку и принимают решение по критерию "хи-квадрат".
2) Проверка появления серий из одинаковых цифр
Обозначим через nL число серий одинаковых подряд цифр длины L. Проверять надо все L от 1 до m, где m - это заданное пользователем число: максимально встречающееся число одинаковых цифр в серии.
В примере "24633899915467766618" обнаружены 2 серии длиной в 2 (33 и 77), то есть n2 = 2 и 2 серии длиной в 3 (999 и 666), то есть n3 = 2.
Вероятность появления серии длиной в L равна: pL = 9 · 10-L (теоретическая). То есть вероятность появления серии длиной в один символ равна: p1 = 0.9 (теоретическая). Вероятность появления серии длиной в два символа равна: p2 = 0.09 (теоретическая). Вероятность появления серии длиной в три символа равна: p3 = 0.009 (теоретическая).
Например, вероятность появления серии длиной в один символ равна pL = 0.9, так как всего может встретиться один символ из 10, а всего символов 9 (ноль не считается). А вероятность того, что подряд встретится два одинаковых символа "XX" равна 0.1 · 0.1 · 9, то есть вероятность 0.1 того, что в первой позиции появится символ "X", умножается на вероятность 0.1 того, что во второй позиции появится такой же символ "X" и умножается на количество таких комбинаций 9.
Частость появления серий подсчитывается по ранее разобранной нами формуле "хи-квадрат" с использованием значений pL.
Примечание: генератор может быть проверен многократно, однако проверки не обладают свойством полноты и не гарантируют, что генератор выдает случайные числа. Например, генератор, выдающий последовательность 12345678912345…, при проверках будет считаться идеальным, что, очевидно, не совсем так.
В заключение отметим, что третья глава книги Дональда Э. Кнута "Искусство программирования" (том 2) полностью посвящена изучению случайных чисел. В ней изучаются различные методы генерирования случайных чисел, статистические критерии случайности, а также преобразование равномерно распределенных случайных чисел в другие типы случайных величин. Изложению этого материала уделено более двухсот страниц.
Лекция 3. Моделирование случайного события
Моделирование полной группы несовместных событий
Моделирование случайного события
Начнем с самого простого. Используем наше умение генерировать случайные числа для имитации выпадения случайных событий.
Случайное событие подразумевает, что у некоторого события есть несколько исходов и то, который из исходов произойдет в очередной раз, определяется только его вероятностью. То есть исход выбирается случайно с учетом его вероятности.
Например, допустим, что нам известна вероятность выпуска бракованных изделий Pб = 0.1. Смоделировать выпадение этого события можно, разыграв равномерно распределенное случайное число из диапазона от 0 до 1 и установив, в какой из двух интервалов (от 0 до 0.1 или от 0.1 до 1) оно попало (см. рис. 3.1). Если число попадает в диапазон (0; 0.1), то выпущен брак, то есть событие произошло, иначе - событие не произошло (выпущено кондиционное изделие). При значительном числе экспериментов частота попадания чисел в интервал от 0 до 0.1 будет приближаться к вероятности P = 0.1, а частота попадания чисел в интервал от 0.1 до 1 будет приближаться к Pк = 0.9.
Рис. 3.1 Схема использования генератора случайных чисел для имитации случайного события
Фрагмент алгоритма представлен на рис. 3.2
Рис. 3.2 Блок-схема алгоритма имитации случайного события
Заметим, что не важно, как вы расположите на отрезке [0; 1] интервал Pб - в начале или в конце, поскольку метод Монте-Карло учитывает только частоту попадания случайных точек в интервал, а она зависит только от величины интервала и не зависит от его месторасположения.
Моделирование полной группы несовместных событий
События называются несовместными, если вероятность появления этих событий одновременно равна 0. Отсюда следует, что суммарная вероятность группы несовместных событий равна 1.
Обозначим через a1, a2, …, an события, а через P1, P2, …, Pn - вероятности появления отдельных событий.
Так как события несовместны, то сумма вероятностей их выпадения равна 1:
P1 + P2 + … + Pn = 1.
Снова используем для имитации выпадения одного из событий генератор случайных чисел, значение которых также всегда находится в диапазоне от 0 до 1. Отложим на единичном интервале [0; 1] отрезки P1, P2, …, Pn. Понятно, что в сумме отрезки составят точно единичный интервал. Точка, соответствующая выпавшему числу из ГСЧ на этом интервале, укажет на один из отрезков. Соответственно в большие отрезки случайные числа будут попадать чаще (вероятность появления этих событий больше!), в меньшие отрезки - реже (см. рис. 3.3).
Рис. 3.3 Схема генерации несовместных случайных событий с помощью генератора случайных чисел
На рис. 3.4 показана блок-схема, которая реализует описанный алгоритм. Алгоритм определяет с помощью фильтра, построенного в виде последовательности условных операций (IF), в какой из интервалов - от 0 до P1, от P1 до (P1 + P2), от (P1 + P2) до (P1 + P2 + P3) и так далее - попало число, сгенерированное генератором случайных чисел. Если число попало в какой-то из интервалов (что произойдет всегда и обязательно), то это соответствует выпадению связанного с ним события.
Рис. 3.4 Блок-схема алгоритма имитации случайных несовместных событий
Пример с возможным исходом четырех несовместных случайных событий.
Промоделируем выпадение последовательности событий - будем выбирать из колоды карт наугад карту (определять ее масть). Карты в колоду возвращать не будем.
В колоде 36 карт четырех мастей по 9 карт каждой масти. Интервал от 0 до 1 разделим на равные четыре части: [0.00-0.25], [0.25-0.50], [0.50-0.75], [0.75-1.00]. Первая часть будет соответствовать картам масти червей (Ч), вторая - картам масти пик (П), третья - картам масти виней (В), четвертая - бубей (Б).
Взять случайное равномерно распределенное число в интервале от 0 до 1 из таблицы случайных чисел или стандартного ГСЧ. Пусть, например, это будет число 0.597. Данное число попадает в третий интервал, соответствующий масти В. Произошло случайное событие: "Масть выпавшей карты - В".
Подобные документы
Распределение дискретной случайной величины по геометрическому закону распределения, проверка теоремы Бернулли на примере моделирования электрической схемы. Математическое моделирование в среде Turbo Pascal. Теоретический расчёт вероятности работы цепи.
контрольная работа [109,2 K], добавлен 31.05.2010Теория вероятностей и математическая статистика являются науками о методах количественного анализа массовых случайных явлений. Множество значений случайной величины называется выборкой, а элементы множества – выборочными значениями случайной величины.
реферат [77,8 K], добавлен 26.12.2008Моделирование случайной величины, распределённой по нормальному закону. Построение доверительных интервалов для математического ожидания и дисперсии, соответствующих доверительной вероятности. Оценка статистических характеристик случайного процесса.
курсовая работа [744,3 K], добавлен 07.06.2010Назначение и принципы действия корреляционно-экстремальной навигационной системы, особенности ее программно-аппаратной реализации, целесообразность статистического моделирования. Описание технологического процесса разработки и отладки программы.
магистерская работа [1,5 M], добавлен 06.12.2013Формулировка теоремы Бернулли, проверка ее с помощью программы. Моделирование случайной величины методом кусочной аппроксимации. График распределения Коши, построение гистограммы и нахождения числовых характеристик, составление статистического ряда.
курсовая работа [226,8 K], добавлен 31.05.2010Теорема Бернулли на примере моделирования электросхемы. Моделирование случайной величины, имеющей закон распределения модуля случайной величины, распределенной по нормальному закону. Проверка критерием Х2: имеет ли данный массив закон распределения.
курсовая работа [2,3 M], добавлен 31.05.2010События и случайные величины. Функция распределения и ее характерные свойства. Сущность и определение основных числовых характеристик случайных величин: математическое ожидание, дисперсия, моменты. Критерии и факторы, влияющие на их формирование.
контрольная работа [118,5 K], добавлен 30.01.2015Понятия теории вероятностей и математической статистики, применение их на практике. Определение случайной величины. Виды и примеры случайных величин. Закон распределения дискретной случайной величины. Законы распределения непрерывной случайной величины.
реферат [174,7 K], добавлен 25.10.2015Построение доверительных интервалов для математического ожидания и дисперсии, соответствующие вероятности. Исследование статистических характеристик случайной величины на основе выбора объема. Теоретическая и эмпирическая плотность распределения.
курсовая работа [594,4 K], добавлен 02.01.2012Особенности выполнения теоремы Бернулли на примере электрической схемы. Моделирование случайной величины по закону распределения Пуассона, заполнение массива. Теория вероятности, понятие ожидания, дисперсии случайной величины и закон распределения.
курсовая работа [29,7 K], добавлен 31.05.2010