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

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

Введение

Случайные процессы - современный и актуальный подход к изучению вероятностных наук. Для более эффективного продвижения в данной области спроектируем, разработаем (реализуем) некоторые концепции из неё на практике.

Таким образом, тема данного дипломного проекта звучит как «Комплекс программ моделирования случайных процессов». Она актуальна, ибо с увеличением системной сложности окружающего нас мира количество (и разнообразие, и скорость их смены) случайных процессов только увеличивается - и нужно уметь с ними работать.

Целью данного труда есть проанализировать теоретические и реализовать практические концепции рассматриваемого подхода. Для её реализации выделим следующие задачи:

.        Повторение (изучение) основ теории случайных процессов.

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

.        Создание программ моделирования случайных последовательностей.

.        Создание программ оценивания моментных функций стационарной случайной последовательности.

.        Создание программы прогнозирования случайной последовательности общего типа.

.        Создание программы, реализующей фильтр Калмана для векторной гауссовской марковской последовательности.

.        Применение разработанных программ в примере прогнозирования случайной последовательности.

Путём напряжённой работы попробуем освоиться в данной области, применяя самые актуальные на сегодняшний день инструменты - превосходную среду математических вычислений Matlab 2016, текстовый редактор Word и другие.

1. Повторение (изучение) основ теории случайных процессов

1.1 Случайный процесс

Случайный процесс (вероятностный процесс, случайная функция, стохастический процесс) в теории вероятностей - семейство случайных величин, индексированных некоторым параметром, чаще всего играющим роль времени или координаты.

Пусть дано вероятностное пространство . Параметризованное семейство  случайных величин , где Т произвольное множество, называется случайной функцией.

Обычно рассматривают два случая:

·        если , то параметр  может интерпретироваться как время. Тогда случайная функция  называется случайным процессом;

·        если множество  дискретно, например , то такой случайный процесс называется случайной последовательностью;

·        если , где , то параметр  может интерпретироваться как точка в пространстве, и тогда случайную функцию называют случайным полем.

Данный подход нестрогий. В частности, термин «случайный процесс» часто используется как безусловный синоним термина «случайная функция».

Классификация:

·        случайный процесс  называется процессом дискретным во времени, если система, в которой он протекает, меняет свои состояния только в моменты времени , число которых конечно или счётно. Случайный процесс называется процессом с непрерывным временем, если переход из состояния в состояние может происходить в любой момент времени;

·        случайный процесс называется процессом с непрерывными состояниями, если значением случайного процесса является непрерывная случайная величина. Случайный процесс называется случайным процессом с дискретными состояниями, если значением случайного процесса является дискретная случайная величина;

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

·        случайная функция называется стационарной в широком смысле, если её математическое ожидание и дисперсия постоянны, а АКФ зависит только от разности моментов времени, для которых взяты ординаты случайной функции. Понятие ввёл А. Я. Хинчин;

·        случайный процесс называется процессом со стационарными приращениями определённого порядка, если вероятностные закономерности такого приращения неизменны во времени. Такие процессы были рассмотрены Ягломом;

·        если ординаты случайной функции подчиняются нормальному закону распределения, то и сама функция называется нормальной;

·        случайные функции, закон распределения ординат которых в будущий момент времени полностью определяется значением ординаты процесса в настоящий момент времени и не зависит от значений ординат процесса в предыдущие моменты времени, называются марковскими;

·        случайный процесс называется процессом с независимыми приращениями, если для любого набора , где , а , случайные величины  независимы в совокупности;

·        если при определении моментных функций стационарного случайного процесса операцию усреднения по статистическому ансамблю можно заменить усреднением по времени, то такой стационарный случайный процесс называется эргодическим;

·        среди случайных процессов выделяют импульсные случайные процессы - случайные процессы, представляющие собой последовательности одиночных импульсов, параметры которых случайно меняются от импульса к импульсу. Под параметрами импульса понимают время его появления, его длительность и параметры, которые определяют его форму.

.2 Статистическое моделирование

Статистическое моделирование - базовый метод моделирования, заключающийся в том, что модель испытывается множеством случайных сигналов с заданной плотностью вероятности. Целью является статистическое определение выходных результатов. В основе статистического моделирования лежит метод Монте-Карло. Отметим, что имитацию используют тогда, когда другие методы применить невозможно.

Рассмотрим метод Монте-Карло на примере вычисления интеграла, значение которого аналитическим способом найти не удается.

Задача 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.1 случаи для N = 1, N = 2, N = 3 - сначала значения частости нельзя назвать достоверными. Попробуем построить график зависимости Pо от N - и посмотрим, как меняется частость выпадения орла в зависимости от количества проведенных опытов. Разумеется, при различных экспериментах будут получаться разные таблицы и, следовательно, разные графики. На рисунке 1.4 показан один из вариантов.

Рисунок 1.4 - Экспериментальная зависимость частости появления случайного события от количества наблюдений и ее стремление к теоретической вероятности

Сделаем некоторые выводы.

.        Видно, что при малых значениях 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); по крайней мере, у них не изменяется первая цифра после запятой, которой мы обязаны доверять по условиям задачи.

.        Причиной такого поведения кривой является действие центральной предельной теоремы. Пока здесь мы сформулируем ее в самом простом варианте «Сумма случайных величин есть величина неслучайная». Мы использовали среднюю величину Pо, которая несет в себе информацию о сумме опытов, и поэтому постепенно эта величина становится все более достоверной.

.        Если проделать еще раз этот опыт сначала, то, конечно, его результатом будет другой вид случайной кривой. И ответ будет другим, хотя примерно таким же. Проведем целую серию таких экспериментов (рисунок 1.5). Такая серия называется ансамблем реализаций. Какому же ответу в итоге следует верить? Ведь они, хоть и являются близкими, все же разнятся. На практике поступают по-разному. Первый вариант - вычислить среднее значение ответов за несколько реализаций (таблица 1.2).

Рисунок 1.5 - Экспериментально снятый ансамбль случайных зависимостей частости появления случайного события от количества наблюдений

Мы поставили несколько экспериментов и определяли каждый раз, сколько необходимо было сделать опытов, то есть Nкрэ. Было проделано 10 экспериментов, результаты которых были сведены в таблице 1.2. По результатам 10-ти экспериментов было вычислено среднее значение Nкрэ.

Таблица 1.2. Экспериментальные данные необходимого количества бросков монеты для достижения точности ε = 0.1 при вычислении вероятности выпадения орла

Таким образом, проведя 10 реализаций разной длины, мы определили, что достаточно в среднем было сделать 1 реализацию длиной в 94 броска монеты.

Еще один важный факт. Внимательно рассмотрим график на рисунке 1.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

Итак, рассмотренные графики подтверждают приведенную выше оценку:

Заметим, что оценок точности может быть несколько.

Пример 1.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. Оценка точности ответа

Поскольку в ответе все еще меняется значение второго разряда, то возможная неточность составляет пока больше 10%. Точность расчета может быть увеличена с ростом числа испытаний (рисунок 1.10).

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

2. Изучение алгоритмов компьютерного моделирования случайных последовательностей

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)

За эталон генератора случайных чисел (ГСЧ) принят такой генератор, который порождает последовательность случайных чисел с равномерным законом распределения в интервале (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. То есть в идеальном случае в каждый интервал попадает одинаковое число точек:

= N/k,

где N - общее число точек, k - количество интервалов, i = 1, …, k.

Рисунок 2.3 - Частотная диаграмма выпадения случайных чисел, порождаемых идеальным генератором теоретически

Следует помнить, что генерация произвольного случайного числа состоит из двух этапов:

·        генерация нормализованного случайного числа (то есть равномерно распределенного от 0 до 1);

·        преобразование нормализованных случайных чисел ri в случайные числа xi, которые распределены по необходимому пользователю (произвольному) закону распределения или в необходимом интервале.

Генераторы случайных чисел по способу получения чисел делятся на:

·        физические;

·        табличные;

·        алгоритмические.

Примером физических ГСЧ могут служить: монета («орел» - 1, «решка» - 0); игральные кости; поделенный на секторы с цифрами барабан со стрелкой; аппаратурный генератор шума (ГШ), в качестве которого используют шумящее тепловое устройство, например, транзистор (рисунки 2.4-2.5).

Рисунок 2.4 - Схема аппаратного метода генерации случайных чисел

Рисунок 2.5 - Диаграмма получения случайных чисел аппаратным методом

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

В персональных компьютерах авторы программных ГСЧ используют гораздо более быстрые источники энтропии, такие, как шум звуковой карты или счётчик тактов процессора. Примеры ГСЧ и источников энтропии отображены в таблице 2.1.

Таблица 2.1. Примеры генераторов случайных чисел и их источников энтропии

Табличные ГСЧ в качестве источника случайных чисел используют специальным образом составленные таблицы, содержащие проверенные некоррелированные, то есть никак не зависящие друг от друга, цифры. В таблице 2.2 приведен небольшой фрагмент такой таблицы. Обходя таблицу слева направо сверху вниз, можно получать равномерно распределенные от 0 до 1 случайные числа с нужным числом знаков после запятой (в нашем примере мы используем для каждого числа по три знака). Так как цифры в таблице не зависят друг от друга, то таблицу можно обходить разными способами, например, сверху вниз, или справа налево, или, скажем, можно выбирать цифры, находящиеся на четных позициях.

Таблица 2.2. Случайные цифры. Равномерно распределенные от 0 до 1 случайные числа

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

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

+ 1 = f(ri).

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

Достоинством данных ГСЧ является быстродействие; генераторы практически не требуют ресурсов памяти, компактны. Недостатки: числа нельзя в полной мере назвать случайными, поскольку между ними имеется зависимость, а также наличие периодов в последовательности квазислучайных чисел.

Важными алгоритмическими методами получения ГСЧ есть:

·        метод серединных квадратов;

·        метод серединных произведений;

·        метод перемешивания;

·        линейный конгруэнтный метод.

Проверка качества работы генератора.

От качества работы ГСЧ зависит качество работы всей системы и точность результатов. Поэтому случайная последовательность, порождаемая ГСЧ, должна удовлетворять целому ряду критериев.

Осуществляемые проверки бывают двух типов:

·        проверки на равномерность распределения;

·        проверки на статистическую независимость.

Проверки на равномерность распределения.

) ГСЧ должен выдавать близкие к следующим значения статистических параметров, характерных для равномерного случайного закона:

- математическое ожидание;

 - дисперсия;

 - среднеквадратичное отклонение.

) Частотный тест.

Частотный тест позволяет выяснить, сколько чисел попало в интервал (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.6).

Рисунок 2.6 - Частотная диаграмма идеального ГСЧ в случае проверки его на частотный тест

Также необходимо учитывать, что количество чисел, попавших в интервал (0; 0,5), должно быть примерно равно количеству чисел, попавших в интервал (0,5; 1).

) Проверка по критерию «хи-квадрат».

Критерий «хи-квадрат» (χ2-критерий) - это один из самых известных статистических критериев; он является основным методом, используемым в сочетании с другими критериями. Критерий «хи-квадрат» был предложен в 1900 году Карлом Пирсоном. Его замечательная работа рассматривается как фундамент современной математической статистики.

Для нашего случая проверка по критерию «хи-квадрат» позволит узнать, насколько созданный нами реальный ГСЧ близок к эталону ГСЧ, то есть удовлетворяет ли он требованию равномерного распределения или нет.

Проверки на статистическую независимость.

1) Проверка на частоту появления цифры в последовательности.

Рассмотрим пример. Случайное число 0,2463389991 состоит из цифр 2463389991, а число 0,5467766618 состоит из цифр 5467766618. Соединяя последовательности цифр, имеем: 24633899915467766618.

Понятно, что теоретическая вероятность pi выпадения i-ой цифры (от 0 до 9) равна 0,1.

Далее следует вычислить частоту появления каждой цифры в выпавшей экспериментальной последовательности. Например, цифра 1 выпала 2 раза из 20, а цифра 6 выпала 5 раз из 20.

Далее считают оценку и принимают решение по критерию «хи-квадрат».

) Проверка появления серий из одинаковых цифр.

Обозначим через 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) полностью посвящена изучению случайных чисел. В ней изучаются различные методы генерирования случайных чисел, статистические критерии случайности, а также преобразование равномерно распределенных случайных чисел в другие типы случайных величин. Изложению этого материала уделено более двухсот страниц.

2.2 Моделирование случайной величины с заданным законом распределения

Большей информативностью, по сравнению с такими статистическими характеристиками как математическое ожидание, дисперсия, для инженера обладает закон распределения вероятности случайной величины X. Представим, что X принимает случайные значения из некоторого диапазона. Например, X - диаметр вытачиваемой детали. Диаметр может отклоняться от запланированного идеального значения под влиянием различных факторов, которые нельзя учесть, поэтому он является случайной слабо предсказуемой величиной. Но в результате длительного наблюдения за выпускаемыми деталями можно отметить, сколько деталей из 1000 имели диаметр X1 (обозначим NX1), сколько деталей имели диаметр X2 (обозначим NX2) и так далее. В итоге можно построить гистограмму частости диаметров, откладывая для X1 величину NX1/1000, для X2 величину NX2/1000 и так далее. (Обратим внимание, если быть точным, NX1 - это число деталей, диаметр которых не просто равен X1, а находится в диапазоне от X1 - Δ/2 до X1 + Δ/2, где Δ = X1 - X2). Важно, что сумма всех частостей будет равна 1 (суммарная площадь гистограммы неизменна). Если X меняется непрерывно, опытов проведено очень много, то в пределе N -> ∞ гистограмма превращается в график распределения вероятности случайной величины. На рисунке 2.7, а показан пример гистограммы дискретного распределения, а на рисунке 2.7, б показан вариант непрерывного распределения случайной величины.

Рисунок 2.7 - Сравнение дискретного и непрерывного законов распределения случайной величины

В нашем примере закон распределения вероятности случайной величины показывает насколько вероятно то или иное значение диаметра выпускаемых деталей. Случайной величиной является диаметр детали.

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

В числе наиболее значимых методов можно назвать:

·        метод ступенчатой аппроксимации;

·        метод усечения;

·        метод взятия обратной функции.

Каждый метод обладает своими положительными и отрицательным сторонами. Среди их реализаций можно выделить достаточно любопытный случай (в большей степени опирающийся на метод ступенчатой аппроксимации).

За неимением места подробно его излагать не будем, перечислим лишь его достоинства:

1.   В качестве функции распределения может быть задана произвольная зависимость f(x), от нее не требуется аналитического задания, непрерывности, дифференцируемости, хорошей интегрируемости (хотя требуется ограниченность на интересующем нас интервале).

2.      Область распределения может быть ограничена любым интервалом.

.        Функция не требует нормировки. Это особенно удобно на ограниченном (по тем или иным соображениям) интервале.

Если же системно (применяя достижения алгебры и численных методов, опирающиеся на твёрдую поверхность кибернетики) подходить к каждому конкретному распределению, то на сегодняшний день можно отметить работы Александра Самарина (The_Freeman) [19, 20], который подвёл итоги применяющимся сегодня (в том числе в САПР вроде Matlab) методам генерации случайных величин с заданным законом распределения. Для примера коротко опишем наиболее часто используемое распределение - нормальное (рисунок 2.8).

Рисунок 2.8 - Функция нормального распределения

Метод инверсии потребует вычисления обратной функции ошибок. Если не использовать специальные аппроксимирующие функции, сложно и невероятно долго. Для нормальной величины существует метод Бокса-Мюллера. Он довольно прост и широко распространен. Его явный недостаток - это вычисление трансцендентных функций. Существует также полярный метод, помогающий избежать подсчета синуса и косинуса. Но мы все еще должны считать логарифм и корень из него. Куда быстрее работает красиво названный метод Ziggurat, придуманный Джорджем Марсалья, автором того же полярного метода.

Полярный метод - это пример выборки с отклонением (acceptance-rejection sampling). Буквально, нужно сгенерировать величину и принять ее, если она подходит, иначе - отклонить и сгенерировать еще раз. Основной пример: генерирация случайной величины с плотностью f(x), однако это слишком сложно сделать простым методом инверсии. Зато, вы можете сгенерировать случайную величину с плотностью g(x), не очень сильно отличающейся от f(x). В таком случае берём наименьшую константу M, такую что M > 1 и почти всюду f(x) < Mg(x), и выполняете следующий алгоритм:

.        Генерируем случайную величину X из распределения g(x) и случайную величину U, равномерно распределенную от 0 до 1.

.        Если U < f(X) / Mg(X) - возвращаем X, иначе - повторяем заново.

Общая проблема выборки с отклонением заключается в подборе такой случайной величины с плотностью распределения g(x), чтобы отклонений было как можно меньше. Для решения этой проблемы существует множество расширений. Сам же метод является основой для почти все последующих алгоритмов, включая Ziggurat. Суть последнего все та же: пытаемся покрыть функцию плотности нормального распределения похожей и более простой функцией и возвращаем величины, попавшие под кривую. Функция своеобразная и напоминает многоступенчатое сооружение, откуда, собственно, и такое название у алгоритма (рисунок 2.9).

Рисунок 2.9 - Схематическая иллюстрация метода Ziggurat

Сам алгоритм можно ускорить различными хитростями работы с 32-битными целыми, убирая ненужные перемножения. Для большей информации можно обратиться к статье George Marsaglia and Wai Wan Tsang «The Ziggurat Method for Generating Random Variables».

Аналогично рассматриваются остальные распределения - непрерывные и дискретные. Имея равномерное и нормальное распределения, часто можно применить формулы преобразования их в требуемое, но, как было отмечено выше, для каждого конкретного случая (в реалиях имеющихся Тьюринг-полных императивных процедурных (наиболее подходящих для быстрых вычислений) языков программирования на ПК с архитектурой фон Неймана) есть свой, особый наиболее оптимальный сегодня алгоритм.

3. Создание программ моделирования случайных последовательностей

.1 Векторная гауссовская марковская последовательность

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

Наиболее известным примером случайного процесса являются флуктуационные (дробовые и тепловые) шумы в радиотехнических устройствах. При наблюдении шумового напряжения на выходах идентичных приборов обнаруживается, что функции, описывающие во времени эти напряжения, различны. Задача теории случайных процессов заключается в отыскании вероятностных закономерностей, которые связывают эти различные функции, описывающие одно и то же физическое явление (например, флуктуационный шум).

Количественно случайный процесс описывается случайной функцией времени ξ(t), которая в любой момент времени t может принимать различные значения с заданным распределением вероятностей. Таким образом, для любого t = ti и значение ξi = ξ(ti) является случайной величиной. Случайный процесс (случайная функция времени) определяется совокупностью функций времени и законами, характеризующими свойства совокупности.

Каждая из функций этой совокупности называется реализацией случайного процесса. Реализация случайной функции ξ(t) обозначается через ξi(t), где параметр i может быть любым действительным числом. Детерминированный процесс имеет одну единственную реализацию, описываемую заданной функцией времени.

Таким образом, случайный процесс - это не одна функция (в классическом понимании), а, в общем случае, бесконечное множество реализаций (ансамбль реализаций). В этом смысле случайная функция есть некая обобщенная функция. Нельзя изобразить случайный процесс одной кривой, пусть даже очень замысловатой. Можно изобразить лишь отдельные его реализации. На рисунке 3.1 в качестве примера изображены реализации случайного процесса и сечение.

В зависимости от того, принадлежат ли возможные значения времени t и реализации ξ(t) дискретному множеству чисел или отрезку (а может быть и всей) действительной оси, различают четыре типа случайных процессов:

.        Случайный процесс общего типа: t и ξ(t) могут принимать любые значения на отрезке (или, быть может, на всей) действительной оси.

.        Дискретный случайный процесс: t непрерывно, а величины ξ(t) дискретны.

.        Случайная последовательность общего типа: t дискретно, а ξ(t) может принимать любые значения на отрезке (или на всей) действительной оси.

.        Дискретная случайная последовательность: t и ξ(t) оба дискретны.

Рисунок 3.1 - Ансамбль реализаций и сечение случайного процесса

Случайный процесс называется гауссовским, если совместная плотность вероятности любой конечной совокупности величин нормальная.

Модель гауссовского случайного процесса широко используется в естествознании и технике. В радиотехнике и связи гауссовский случайный процесс является адекватной математической моделью активных и пассивных помех, атмосферных и космических шумов, каналов с замиранием, с многолучевым распространением, групповых сигналов в многоканальных системах. Флуктуационные шумы приемных устройств, обусловленные дробовым эффектом и тепловым движением электронов, также подчиняются нормальному закону распределения. Адекватность модели гауссовского случайного процесса многим реальным помехам и сигналам и ее универсальность объясняются во многих случаях действием центральной предельной теоремы теории вероятностей.

Ковариационная матрица (или матрица ковариаций) в теории вероятностей - это матрица, составленная из попарных ковариаций элементов одного или двух случайных векторов.

Ковариационная матрица случайного вектора - квадратная симметрическая неотрицательно определенная матрица, на диагонали которой располагаются дисперсии компонент вектора, а внедиагональные элементы - ковариации между компонентами.

Ковариационная матрица случайного вектора является многомерным аналогом дисперсии случайной величины для случайных векторов. Матрица ковариаций двух случайных векторов - многомерный аналог ковариации между двумя случайными величинами.

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

Взаимная ковариационная функция случайных процессов - это функция двух переменных t и u, равная математическому ожиданию произведения случайных процессов, взятых в любые моменты времени t и u из областей определения этих случайных процессов.

Разведаем доступные средства новейшего IDE Matlab 2016. Среди них для наших целей может быть использована в числе прочих функция inv(), которая обращает матрицу. «Вокруг неё» и формул теории вероятности и построим нашу программу.

Если ещё конкретнее, то последовательность необходимых действий будет следующей:

1.      Устанавливаются начальные данные:

·        интервал дискретизации;

·        длительность;

·        функция корреляции (произвольная);

·        число отсчетов.

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

.        После генерации четырех первых значений генерируется матрица связи для генерации последующих.

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

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

Содержимое файла GaussVectorSequense.m

function GaussVectorSequense(filename,a,b,N)

if nargin < 4 %по умолчанию

filename='gaussVector'

N=1000;

a=2;

b=4/2*pi;

end

del=1/2; % интервал дискретизации

T=2.5; % длительность реализаций

t=0:del:T;

r=2*exp(-t)-exp(-2*t); % функция корреляции

n=length(t); % число отсчётов

for i=1:n

for j=1:n

b(i,j)=r(abs(i-j)+1);

end

end

BB=b ; % исходная корреляционная матрица

D=inv(BB) ;% матрица точности

D(1,6)=0;

D(6,1)=0;

B=inv(D) ; % аппроксимирующая корреляционная матрица

% корреляционные матрицы первых значений

B3=[B(1,1) B(1,2)

B(2,1) B(2,2)];

B4=[B(1,1) B(1,2) B(1,3)

B(2,1) B(2,2) B(2,3)

B(3,1) B(3,2) B(3,3)];

B5=[B(1,1) B(1,2) B(1,3) B(1,4)

B(2,1) B(2,2) B(2,3) B(2,4)

B(3,1) B(3,2) B(3,3) B(3,4)

B(4,1) B(4,2) B(4,3) B(4,4)];

x(1)=randn*sqrt(B(1,1)) ; % первое значение последовательности

c=B(2,1)/B(1,1);

si=sqrt(B(2,2)-c*B(1,2));

x(2)=c*x(1)+randn*si ; % второе значение последовательности

C=[B(3,1) B(3,2)]*inv(B3);

si=sqrt(B(3,3)-C*[B(1,3);B(2,3)]);

x(3)=C(1)*x(1)+C(2)*x(2)+randn*si ; % третье значение последовательности

C=[B(4,1) B(4,2) B(4,3)]*inv(B4);

si=sqrt(B(4,4)-C*[B(1,4);B(2,4);B(3,4)]);

x(4)=C(1)*x(1)+C(2)*x(2)+C(3)*x(3)+randn*si; % четвёртое значение

C=[B(5,1) B(5,2) B(5,3) B(5,4)]*inv(B5) ;% матрица связи

s=sqrt(B(5,5)-C*[B(1,5);B(2,5);B(3,5);B(4,5)]); % ср.квадр.отклон.

X=[x(1);x(2);x(3);x(4)];

fid = fopen(filename, 'w');

if fid == -1

error('File '+filename+'is not opened');

end

fprintf(fid,'%f ',x(1));

fprintf(fid,'%f ',x(2));

fprintf(fid,'%f ',x(3));

fprintf(fid,'%f ',x(4));

sw = ['' sprintf('\n')];

for i=5:N

x(i)=C*X+s*randn; % следующее значение

fprintf(fid,'%f ',x(i));

if mod(i,10)==0

fprintf(fid,'%s ',sw);

end

for j=1:4

X(j)=x(i-4+j); % сдвиг предыдущих четырёх

end

end

plot(x); % марковская последовательность

Содержимое файла GaussVectorSequense_caller.m

answer=GaussVectorSequense('GaussVector',2,4/2*pi,1000) ;(answer);

Результаты работы отобразим на рисунках 3.2 и 3.3.

Рисунок 3.2 - Текстовая часть результатов работы программы

Рисунок 3.3 - Графическая часть результатов работы программы

.2 Скалярная гауссовская последовательность общего типа (с произвольной ковариационной функцией)

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

Пример: броуновское движение. Частицы (молекулы) газа или жидкости в отсутствие внешних влияний находятся в постоянном, хаотическом движении, интенсивность которого зависит только от температуры и плотности. Частица в случайные моменты времени сталкивается с молекулами окружающей среды и меняет при этом свою скорость и направление.

Будем следить, как изменяется с течением времени одна из координат избранной частицы, допустим горизонтальная координата. При этом можно не учитывать силу тяжести; на частицу действуют только систематическая сила трения об окружающую среду и случайная сила толчков. Считая компоненту силы трения пропорциональной х-й компоненте скорости и пренебрегая силой инерции, получим следующее уравнение движения частицы:

γ = n(t)

где γ - коэффициент трения; n(t) - сила случайных толчков вдоль оси х.

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

Смещения частицы для двух не перекрывающихся интервалов времени, значительно превышающих время свободного пробега, независимы. Если взять три момента времени (рисунок 3.4).

> t’ > t1,

причем промежутки много больше среднего времени свободного пробега, то поведение частицы на интервале (t’, t2) не будет зависеть от того, что происходило с частицей до момента t’. Это характерное свойство Марковских процессов. Такой случайный процесс можно назвать процессом без последействия.

Рисунок 3.4 - Траектории броуновской частицы

Допустим, что нам точно известна координата х’ в момент t’. Вследствие случайного характера воздействующей силы n(t), возможные значения координаты х2 в момент времени t2 различны и образуют некоторый «ансамбль». Мы можем говорить об условной вероятности F(х2, t2 | х’, t’) того, что если в момент времени t’ координата равна х’, то в момент времени t2 частица будет иметь координату, заключенную в промежутке (х2, х2 + dx2).

Условная вероятность F(х2, t2 | х’, t’) характеризует вероятность перехода частицы из состояния х’ в состояние х2 за время между t’ и t2 и называется вероятностью перехода. Если частица в момент времени t’ может иметь различные значения координаты х’ с вероятностью F(х’, t’), то двумерная плотность вероятности равна:

(х’, х2; t’, t2) = f(х’, t’) f(х2, t2 | х’, t’).

Нормальные процессы и марковские процессы в общем случае частично «перекрываются».

Воспользуемся встроенными возможностями новейшего математического («матричного») САПР Matlab 2016. Среди них для наших целей выделяется функция normrnd(), которая возвращает элементы нормального распределения. «Вокруг неё» и построим программу.

Если ещё конкретнее, то последовательность необходимых действий будет следующей:

1.      Пользователь указывает имя выходного файла, МО, дисперсию и число элементов для генерации

.        С помощью функции randn() генерируется нормально распределенная завтравка для генератора, и с помощью классической формулы (плотности) для нормального распределения генерируются выходные элементы, которые записываются в файл.

.        Далее последовательно отрисовывается на графике.

Содержимое файла GaussScalarSequense.mGaussScalarSequense(filename,a,sigma,N)

fid = fopen(filename, 'w'); % открытие файла на запись

if fid == -1 % проверка корректности открытия

error('File '+filename+'is not opened');

end

if nargin < 4 %по умолчанию

N=100;

a=0;

sigma=1;

end

% x=1:1:N;

arr=[0,0];%набор значений случайной величины

i=1;

sw = ['' sprintf('\n')];

while i < N+1

rng('shuffle');

rnd=randn;

arr(i)=((1/sqrt(2*pi*sigma2))*exp(-((rnd-a)2/(2*sigma2))));

fprintf(fid,'%f ',arr(i));

i=i+1;

if mod(i,10)==0

fprintf(fid,'%s ',sw);

end

end

fclose(fid);

i=1:1:N;

plot(i,arr(i))

end

Содержимое файла GaussScalarSequense_caller.m

answer=GaussScalarSequense('Gauss1',0,1,100) ;(answer);

Результаты работы отобразим на рисунках 3.6 и 3.7.

Рисунок 3.6 - Текстовая часть результатов работы программы

Рисунок 3.7 - Графическая часть результатов работы программы

4. Создание программ оценивания моментных функций стационарной случайной последовательности

4.1 Оценивание математического ожидания

Математическое ожидание случайной величины X (обозначается M(X) или реже E(X)) характеризует среднее значение случайной величины (дискретной или непрерывной). Математическое ожидание - это первый начальный момент заданной СВ.

Математическое ожидание относят к так называемым характеристикам положения распределения (к которым также принадлежат мода и медиана). Эта характеристика описывает некое усредненное положение случайной величины на числовой оси. Скажем, если матожидание случайной величины - срока службы лампы, равно 100 часов, то считается, что значения срока службы сосредоточены (с обеих сторон) от этого значения (с тем или иным разбросом, о котором уже говорит дисперсия).

Математическое ожидание дискретной случайной величины Х вычисляется как сумма произведений значений xi, которые принимает СВ Х, на соответствующие вероятности pi:

Для непрерывной случайной величины (заданной плотностью вероятностей f(x)), формула вычисления математического ожидания Х выглядит следующим образом:

Воспользуемся встроенными возможностями мощного САПР Matlab 2016. Среди них для наших целей выделяется функция sum(), которая возвращает сумму чисел. Сделаем её нашей отправной точкой, и придём к успеху (путём деления суммы чисел на их количество (подходящий тип оценки матожидания)).

Если ещё конкретнее, то последовательность необходимых действий будет следующей:

.        Для вычисления оценки математического ожидания необходимо считать последовательность и число ее элементов.

.        Далее простым суммированием элементов и делением суммы на число элементов вычислить оценку математического ожидания.

Содержимое файла GetMiddleOfGauss.m

function a= GetMiddleOfGauss(fileName,N)

fid = fopen(fileName, 'r'); % открытие файла на запись

if fid == -1 % проверка корректности открытия

error('File'+fileName+' is not opened');

end

if nargin < 2 %по умолчанию

N=100;

end

values =fscanf(fid,'%f');

summ=sum(values);

fclose(fid);

a=summ/N;

Содержимое файла GetMiddleOfGauss_caller.m('Gauss1',0,1,1000) ;=GetMiddleOfGauss('Gauss1',1000) ;('Middle of gausse sequense:');(answer);

Результаты работы отобразим на рисунке 4.1.

Рисунок 4.1 - Результаты работы программы

.2 Оценивание ковариационной функции

Автокорреляция - статистическая взаимосвязь между последовательностями величин одного ряда, взятыми со сдвигом, например, для случайного процесса - со сдвигом по времени.

Данное понятие широко используется в эконометрике. Наличие автокорреляции случайных ошибок регрессионной модели приводит к ухудшению качества МНК-оценок параметров регрессии, а также к завышению тестовых статистик, по которым проверяется качество модели (то есть создается искусственное улучшение качества модели относительно её действительного уровня точности). Поэтому тестирование автокорреляции случайных ошибок является необходимой процедурой построения регрессионной модели.

Коэффициенты автокорреляции также имеют самостоятельное важное значение для моделей временных рядов ARMA.

Чаще всего тестируется наличие в случайных ошибках авторегрессионного процесса первого порядка. Для тестирования нулевой гипотезы, о равенстве коэффициента автокорреляции нулю чаще всего применяют критерий Дарбина-Уотсона. При наличии лаговой зависимой переменной в модели данный критерий неприменим, можно использовать асимптотический h-тест Дарбина. Оба эти теста предназначены для проверки автокорреляции случайных ошибок первого порядка. Для тестирования автокорреляции случайных ошибок большего порядка можно использовать более универсальный асимптотический LM-тест Бройша-Годфри. В данном тесте случайные ошибки не обязательно должны быть нормально распределены. Тест применим также и в авторегрессионных моделях (в отличие от критерия Дарбина-Уотсона).

Для тестирования совместной гипотезы о равенстве нулю всех коэффициентов автокорреляции до некоторого порядка можно использовать Q-тест Бокса-Пирса или Q-тест Льюнга-Бокса

Автокорреляционная функция показывает зависимость автокорреляции от величины сдвига во времени. При этом предполагается стационарность временного ряда, означающая в том числе независимость автокорреляций от момента времени. Анализ автокорреляционной функции (вместе с частной автокорреляционной функцией) позволяет проводить идентификацию порядка ARMA-моделей.

Автоковариационная функция - произведение автокорреляционной функции на дисперсию. Обычно используется именно автокорреляционная функция, так как она безразмерна, то есть независима от масштаба изменения временных рядов.

Воспользуемся встроенными возможностями современного математического («матричного») САПР Matlab 2016. Среди них для наших целей выделяется функция cov(), которая возвращает ковариационную функцию двух распределений. Продвинемся вперёд, опираясь на неё.

Если ещё конкретнее, то последовательность необходимых действий будет следующей:

.        Для вычисления ковариации двух последовательностей необходимо их считать из файла и подсчитать сумму элементов для каждой последовательности.

.        Затем подсчитать оценку МО для каждой последовательности.

.        После чего, воспользовавшись формулой оценивания ковариации вычислить ее.

Содержимое файла GetGaussCovariation.mcovar= GetGaussCovariation(fileName1,fileName2,N)

fid1 = fopen(fileName1, 'r'); % открытие файла на чтение

fid2 = fopen(fileName2, 'r');

if fid1 == -1 % проверка корректности открытия

error('File'+fileName1+' is not opened');

end

if fid2 == -1 % проверка корректности открытия

error('File'+fileName2+' is not opened');

end

if nargin < 3 %по умолчанию

N=100;

end

values1 =fscanf(fid1,'%f');

values2 =fscanf(fid2,'%f');

fclose(fid1);

fclose(fid2);

summ1=sum(values1);

summ2=sum(values2);

a1=summ1/N;

a2=summ2/N;

i=1;

total=0;

while i<N+1

total=total+(values1(i)-a1)*(values2(i)-a2);

i=i+1;

end

covar=total/N;

Содержимое файла GetGaussCovariation_caller.m('Gauss1',0,1,100) ;('Gauss2',0,2,100) ;=GetGaussCovariation('Gauss1','Gauss2',100) ;('Covariation of gausse sequenseS:');(answer);

Результаты работы отобразим на рисунках 4.2 и 4.3.

Рисунок 4.2 - Текстовая часть результатов работы программы

Рисунок 4.3 - Графическая часть результатов работы программы

5. Создание программы прогнозирования случайной последовательности общего типа

.1 Короткий теоретический базис

Функция правдоподобия в математической статистике - это совместное распределение выборки из параметрического распределения, рассматриваемое как функция параметра. При этом используется совместная функция плотности (в случае выборки из непрерывного распределения) либо совместная вероятность (в случае выборки из дискретного распределения), вычисленные для данных выборочных значений.

Понятия вероятности и правдоподобия тесно связаны. Сравните два предложения:

«Какова вероятность выпадения 12 очков в каждом из ста бросков двух костей?»

«Насколько правдоподобно, что кости не шулерские, если из ста бросков в каждом выпало 12 очков?»

Если распределение вероятности зависит от параметра, то с одной стороны можно рассматривать вероятность некоторых событий при заданном параметре, а с другой стороны - вероятность заданного события при различных значениях параметра. То есть в первом случае имеем функцию, зависящую от события, а во втором - от параметра при фиксированном событии. Последний вариант является функцией правдоподобия и показывает, насколько правдоподобен выбранный параметр при заданном событии.

Неформально: если вероятность позволяет нам предсказывать неизвестные результаты, основанные на известных параметрах, то правдоподобие позволяет нам оценивать неизвестные параметры, основанные на известных результатах.

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

Метод максимального правдоподобия или метод наибольшего правдоподобия (ММП, ML, MLE - англ. maximum likelihood estimation) в математической статистике - это метод оценивания неизвестного параметра путём максимизации функции правдоподобия. Основан на предположении о том, что вся информация о статистической выборке содержится в функции правдоподобия. Метод максимального правдоподобия был проанализирован, рекомендован и значительно популяризирован Р. Фишером между 1912 и 1922 годами (хотя ранее он был использован Гауссом, Лапласом и другими).

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

Метод максимального правдоподобия соответствует многим известным методам оценки в области статистики. Например, вы интересуетесь таким антропометрическим параметром, как рост жителей России. Предположим, у вас имеются данные о росте некоторого количества людей, а не всего населения. Кроме того, предполагается, что рост является нормально распределённой величиной с неизвестной дисперсией и средним значением. Среднее значение и дисперсия роста в выборке являются максимально правдоподобными к среднему значению и дисперсии всего населения.

Для фиксированного набора данных и базовой вероятностной модели, используя метод максимального правдоподобия, мы получим значения параметров модели, которые делают данные «более близкими» к реальным. Оценка максимального правдоподобия даёт уникальный и простой способ определить решения в случае нормального распределения.

Метод оценки максимального правдоподобия применяется для широкого круга статистических моделей, в том числе:

·        линейные модели и обобщённые линейные модели;

·        факторный анализ;

·        моделирование структурных уравнений;

·        многие ситуации, в рамках проверки гипотезы и доверительного интервала формирования;

·        дискретные модели выбора.

.2 Проектирование

Программа прогнозирования стоимости на электричество, файл prices.mat содержит порядка 50000 значений: дата-час-стоимость.

Прогнозируется на 24 шага вперед, иными словами на сутки.

В прогнозировании используется метод выборки максимального правдоподобия.

Модель схожа с моделью исторического развития по спирали: исторические этапы повторяются, но свойства и характеристики этапов изменяются.

Если рассматривать последовательности временных отметок, то если значений последовательности достаточно много, то практически всегда можно выделить некоторый интервал, который похож на тот, что происходит накануне прогноза.

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

Метод подобия выборок применимы в краткосрочных прогнозах (финансы, энергетика, ...).

Исходные данные в файле prices.mat (открывать в Matlab)

Задействуем средства мощного САПР Matlab 2016. Из них для наших целей можно особо выделитть класс timeseries(), который содержит данные - вектор чисел (часто через равные интервалы) и методы - для идентификации и моделирования закономерностей, а также прогнозирования на их основе данных. Сделаем его нашей отправной точкой, и получим результат.

.3 Разработка

Содержимое файла task5.m

clear; clc;

% <Date> - <Hour> - <Value>prices PRICES_EUR; % Загрузка файла prices.mat= PRICES_EUR;

% Задача прогнозирования= datenum('01.09.2012 23:00:00', 'dd.mm.yyyy hh:MM:ss'); % Момент прогноза = 24; % Горизонт прогнозирования (Forecast)

% Параметры максимального подобия

M = 144;

Step = 24;

% Определяется выбока новой истории

Index = find(TimeSeries(:,1) == datenum(year(T), month(T), day(T)) & TimeSeries(:,2) == hour(T)); length(Index) > 1

fprintf(['Ошибка временного ряда: отметка времени T найдена во временном ряде более 1 раза \n']);

elseif isempty(Index)

fprintf(['Ошибка временного ряда: отметка времени T во временном ряде не найдена \n']);

else

HistNewData = TimeSeries([Index-M+1:Index],:); % Выборка новой истории (HistNewData)

Index = Index - Step * 2;

end

%Определить значения подобия при помощи перебора

k = 1;Index > 2 * M

HistOldData = TimeSeries([Index-M+1 : Index],:);

Likeness(k,1)= Index;

CheckOld = find(HistOldData(:,3) > 0);

CheckNew = find(HistNewData(:,3) > 0);

if isempty(CheckOld) || isempty(CheckNew)

Likeness(k,2) = 0;

else

Likeness(k,2) = abs(corr(HistOldData(:,3), HistNewData(:,3), 'type', 'Pearson'));

end

k = k + 1;

Index = Index - Step; % Возврат по времени назад на шаг Step

end

% 2.2.3. Определить максимум подобия MaxLikeness

MaxLikeness = max(abs(Likeness(:,2)));

IndexLikeness = find(Likeness(:,2) == MaxLikeness);= Likeness(IndexLikeness(1),1);

(['Максимальное подобие = ', num2str(MaxLikeness),'\n'])

% Определить выборку максимального подобия

MSPData = TimeSeries([MSP-M+1 : MSP],:); % Выборка максимального подобия

% Определить базовую выборку (BaseHistData)

HistBaseData = TimeSeries([MSP+1:MSP+P],:);% Базовая выборка

% Коэффициенты alpha1 и alpha0

% Делаем аппроксимацию HistNewData при помощи MSPData по методу наименьших квадратов.

% В данном случае решение находится через обратную матрицу.

X = MSPData(:,3);

X(:,length(X(1,:))+1) = 1; % Добавляем столбец с единичным вектором I

Y = HistNewData(:,3);= X(:,2);= X'* X;= X'* Y;= inv(Xn); = invX * Yn; % Искомые коэффициенты alpha1 и alpha0 являются значениями матрица A

% Собственно, прогнозирование

X = HistBaseData(:,3);

X(:,length(X(1,:))+1) = 1;

Forecast = X * A; % Выборка из 24 прогнозных значения

disp('Спрогнозированный результат:');

disp(Forecast)

5.4 Результаты

Результаты работы отобразим на рисунке 5.1.

Рисунок 5.1 - Результаты работы программы

6. Создание программы, реализующей фильтр Калмана для векторной гауссовской марковской последовательности

6.1 Короткий теоретический базис

Фильтр Калмана - эффективный рекурсивный фильтр, оценивающий вектор состояния динамической системы, используя ряд неполных и зашумленных измерений. Назван в честь Рудольфа Калмана.

Фильтр Калмана широко используется в инженерных и эконометрических приложениях: от радаров и систем технического зрения до оценок параметров макроэкономических моделей. Калмановская фильтрация является важной частью теории управления, играет большую роль в создании систем управления. Совместно с линейно-квадратичным регулятором фильтр Калмана позволяет решить задачу линейно-квадратичного гауссовского управления. Фильтр Калмана и линейно-квадратичный регулятор - возможное решение большинства фундаментальных задач в теории управления.

В большинстве приложений размерность вектора состояния объекта превосходит размерность вектора данных наблюдения. И при этом фильтр Калмана позволяет оценивать полное внутреннее состояние объекта.

Фильтр Калмана предназначен для рекурсивного дооценивания вектора состояния априорно известной динамической системы, то есть для расчёта текущего состояния системы необходимо знать текущее измерение, а также предыдущее состояние самого фильтра. Таким образом, фильтр Калмана, подобно другим рекурсивным фильтрам, реализован во временно́м, а не в частотном представлении, но в отличие от других подобных фильтров, фильтр Калмана оперирует не только оценками состояния, а ещё и оценками неопределенности (плотности распределения) вектора состояния, опираясь на формулу Байеса условной вероятности.

Алгоритм работает в два этапа. На этапе прогнозирования фильтр Калмана экстраполирует значения переменных состояния, а также их неопределенности. На втором этапе, по данным измерения (полученного с некоторой погрешностью), результат экстраполяции уточняется. Благодаря пошаговой природе алгоритма, он может в реальном времени отслеживать состояние объекта (без заглядывания вперед, используя только текущие замеры и информацию о предыдущем состоянии и его неопределенности).

Бытует ошибочное мнение, что для правильной работы фильтра Калмана якобы требуется гауссовское распределение входных данных. В исходной работе Калмана результаты о минимуме ковариации фильтра были получены на базе ортогональных проекций, без предположений о гауссовости ошибок измерений. Затем просто было показано, что для специального случая распределения ошибок по Гауссу фильтр дает точную оценку условной вероятности распределения состояния системы.

Наглядный пример возможностей фильтра - получение оптимальных, непрерывно обновляемых оценок положения и скорости некоторого объекта по результатам временного ряда неточных измерений его местоположения. Например, в радиолокации стоит задача сопровождения цели, определения её местоположения, скорости и ускорения, при этом результаты измерений поступают постепенно и сильно зашумлены. Фильтр Калмана использует вероятностную модель динамики цели, задающую тип вероятного движения объекта, что позволяет снизить воздействие шума и получить хорошие оценки положения объекта в настоящий, будущий или прошедший момент времени.

Существует несколько разновидностей фильтра Калмана, отличающихся приближениями и ухищрениями, которые приходится применять для сведения фильтра к описанному виду и уменьшения его размерности:

·        расширенный фильтр Калмана (EKF, Extended Kalman filter). Сведение нелинейных моделей наблюдений и формирующего процесса с помощью линеаризации посредством разложения в ряд Тейлора;

·        сигма-точечный фильтр Калмана (UKF, Unscented Kalman filter). Используется в задачах, в которых простая линеаризация приводит к уничтожению полезных связей между компонентами вектора состояния. В этом случае «линеаризация» основана на сигма-точечном преобразовании;

·        Ensemble Kalman filter (EnKF). Используется для уменьшения размерности задачи;

·        возможны варианты с нелинейным дополнительным фильтром, позволяющим привести негауссовые наблюдения к нормальным;

·        возможны варианты с «обеляющим» фильтром, позволяющим работать с «цветными» шумами и т. д.

Кроме того, имеются аналоги фильтра Калмана, использующие полностью или частично модель непрерывного времени:

·        фильтр Калмана-Бьюси, в котором и эволюция системы, и измерения имеют вид функций от непрерывного времени;

·        гибридный фильтр Калмана, использующий непрерывное время при описании эволюции системы, и дискретные моменты времени для измерений.

.2 Проектирование

Безусловно самый сложный алгоритм из представленных. Для его реализации нужно:

.        Задать входные данные.

·        K - количество реализаций процесса моделирования фильтра;

·        m - размерность вектора состояния;

·        r - размерность вектора возмущения;

·        n - размерность вектора наблюдений;

·        ro - коэффициент корреляции;

·        dx - дисперсия последовательности;

·        Ft, Gt - системные матрицы уравнения состояний;

·        ht - матрицы связи состояний и наблюдений;

·        Qt - дисперсия шума возмущения;

·        Rt - дисперсия шума наблюдения.

.        После задания входных данных начать цикл моделирования наблюдений.

.        В процессе моделирования наблюдений осуществить проведение оценки с помощью фильтра Калмана в пакетном режиме.

.        Полученные оценки записать в файл.

У нас в распоряжении есть встроенные возможности развитого САПР Matlab 2016. Среди них для наших целей нельзя особо выделить ни одну функцию - для реализации поставленной задачи приходится все формулы набирать вручную. Аналогично предыдущим примерам применим модульную структуру программы, где KalmanFilter_caller.m её запускает, KalmanFilter.m проводит основные операции, а kaldt.m реализует «внутреннюю логику» фильтра.

6.3 Разработка

Содержимое файла kaldt.m

function [xt,xt_,Pt,Pt_]= kaldt(zt,Ft,Gt,ht,Qt,Rt,x10,P10)

[n,t]=size(zt); [m,mf]=size(Ft); [mg,r]=size(Gt); [rq1,rq2]=size(Qt);

[nh,mh]=size(ht); [nr1,nr2]=size(Rt); mx=length(x10);

[mp1,mp2]=size(P10);

xt=zeros(m,t);Pt=zeros(m,m,t);xt_=zeros(m,t+1);

xt_(:,1)=x10; Pt_(:,:,1)=P10;

for cnt=1:t

Vt=Pt_(:,:,cnt)*ht';Ut=ht*Vt+Rt; Ut_=inv(Ut);Wt=Vt*Ut_;

Pt(:,:,cnt)=Pt(:,:,cnt)-Wt*Vt';

xpr=xt_(:,cnt)+Wt*(zt(:,cnt)-ht*xt(:,cnt));xt(:,cnt)=xpr;

xt_(:,cnt+1)=Ft*xpr;

Pt_(:,:,cnt+1)=Ft*Pt(:,:,cnt)*Ft'+Gt*Qt*Gt';

end;

end

Содержимое файла KalmanFilter.m

function KalmanFilter(filenameIn,filenameOut,N)

if nargin < 1 %по умолчанию

filenameIn='gaussVector';

filenameOut='kalmanOut';

N=20;

end

fIn = fopen(filenameIn, 'r'); % открытие файла на запись

if fIn == -1 % проверка корректности открытия

error('File'+filenameIn+' is not opened');

end

values =fscanf(fIn,'%f');

fclose(fIn);

K=length(values);

m=1;r=1;n=1;

%t=1:N;

ro=0.99;dx=1;

Ft=ro;

Gt=sqrt(dx*(1-ro2));

ht=1;

Qt=1;Rt=1;

% Pte=zeros(m,m,N);Pt_e=zeros(m,m,N);

for k=1:K

x=zeros(m,N);values=zeros(n,N);

ut=randn(r,N);

vt=randn(n,N);

x(:,1)=sqrt(dx)*randn;

values(:,1)=ht*x(:,1)+vt(:,1);

for i=1:N-1

x(:,i+1)=Ft*x(:,i)+Gt*ut(:,i);

values(:,i+1)=ht*x(:,i+1)+vt(:,i+1);

end;

x10=0;P10=dx;

[xt,xt_,Pt,Pt_]=kaldt(values,Ft,Gt,ht,Qt,Rt,x10,P10);

end;

disp(xt);

fOut = fopen(filenameOut, 'w');

if fOut == -1

error('File '+filenameOut+'is not opened');

end

sw = ['' sprintf('\n')];

for i=1:length(xt)

fprintf(fOut,'%f ',xt(i));

if mod(i,10)==0

fprintf(fOut,'%s ',sw);

end

end

fclose(fOut);

end

Содержимое файла KalmanFilter_caller.m

answer=KalmanFilter('gaussVector','kalmanFilter',100) ;('Гауссова последователность, подвергнутая фильтру Калмана:');

disp(answer);

%Параметры

%входной файл, выходной файл, количество дискретных отсчетов

.4 Результаты

Результаты работы отобразим на рисунке 6.1.

Рисунок 6.1 - Результаты работы программы

7 Применение разработанных программ в примере прогнозирования случайной последовательности

.1 Проектирование

Применим разработанные программ в примере прогнозирования случайной последовательности. Это несколько затруднительно ввиду их разнородности, но наибольший смысл, как нам представляется, несёт в себе следующая последовательность действий:

1.      Определение выборки максимального подобия MSPData.

.        Определение базовой выборки BaseHistData.

.        Аппроксимация HistNewData при помощи MSPData по методу наименьших квадратов.

.        Вычисление среднего от оригинальной последовательности и от спрогнозированной.

.        Применение фильтра Калмана к оригинальной последовательности и к спрогнозированной.

.        Вычисление ковариации от оригинальной последовательности и спрогнозированной.

7.2 Разработка

Содержимое файла task7_caller.m.m

clear; clc;

% <Date> - <Hour> - <Value>prices PRICES_EUR; % Загрузка файла prices.mat= PRICES_EUR;

% Задача прогнозирования= datenum('01.09.2012 23:00:00', 'dd.mm.yyyy hh:MM:ss'); % Момент прогноза = 24; % Горизонт прогнозирования (Forecast)

% Параметры максимального подобия

M = 144;

Step = 24;

% Определяется выбока новой истории

Index = find(TimeSeries(:,1) == datenum(year(T), month(T), day(T)) & TimeSeries(:,2) == hour(T)); length(Index) > 1

fprintf(['Ошибка временного ряда: отметка времени T найдена во временном ряде более 1 раза \n']);

elseif isempty(Index)

fprintf(['Ошибка временного ряда: отметка времени T во временном ряде не найдена \n']);

else

HistNewData = TimeSeries([Index-M+1:Index],:); % Выборка новой истории (HistNewData)

Index = Index - Step * 2;

end

%Определить значения подобия при помощи перебора

k = 1;Index > 2 * M

HistOldData = TimeSeries([Index-M+1 : Index],:);

Likeness(k,1)= Index;

CheckOld = find(HistOldData(:,3) > 0);

CheckNew = find(HistNewData(:,3) > 0);

if isempty(CheckOld) || isempty(CheckNew)

Likeness(k,2) = 0;

else

Likeness(k,2) = abs(corr(HistOldData(:,3), HistNewData(:,3), 'type', 'Pearson'));

end

k = k + 1;

Index = Index - Step; % Возврат по времени назад на шаг Step

end

% 2.2.3. Определить максимум подобия MaxLikeness

MaxLikeness = max(abs(Likeness(:,2)));

IndexLikeness = find(Likeness(:,2) == MaxLikeness);= Likeness(IndexLikeness(1),1);

% Определить выборку максимального подобия

MSPData = TimeSeries([MSP-M+1 : MSP],:); % Выборка максимального подобия

% Определить базовую выборку (BaseHistData)

HistBaseData = TimeSeries([MSP+1:MSP+P],:);% Базовая выборка

% Коэффициенты alpha1 и alpha0

% Делаем аппроксимацию HistNewData при помощи MSPData по методу наименьших квадратов.

% В данном случае решение находится через обратную матрицу.

X = MSPData(:,3);

X(:,length(X(1,:))+1) = 1; % Добавляем столбец с единичным вектором I

Y = HistNewData(:,3);= X(:,2);= X'* X;= X'* Y;= inv(Xn); = invX * Yn; % Искомые коэффициенты alpha1 и alpha0 являются значениями матрица A

% Собственно, прогнозирование

X = HistBaseData(:,3);

X(:,length(X(1,:))+1) = 1;

Forecast = X * A; % Выборка из 24 прогнозных значения

disp('Спрогнозированная последовательность:')

disp(Forecast)

% 1) Определяем фактические значения= find(TimeSeries(:,1) == datenum(year(T), month(T), day(T)) & TimeSeries(:,2) == hour(T));= TimeSeries([Index : Index+P-1],3);('Реальная последовательность:')

disp(Fact)

filename_origin='original_sequence';_predict='predict_sequence';_origin = fopen(filename_origin, 'w'); _predict = fopen(filename_predict, 'w'); fid_origin == -1

error('File '+filename_origin+'is not opened'); fid_predict == -1

error('File '+filename_predict+'is not opened'); i=1:24

fprintf(fid_origin,'%f ',Fact(i));

fprintf(fid_predict,'%f ',Forecast(i));

end

%Вычисление среднего от оригинальной последовательности и от

%спрогнозированной_origin=GetMiddleOfGauss(filename_origin,24) ;_predict=GetMiddleOfGauss(filename_predict,24) ;('МО реальной последовательности:');

disp(a_origin);

disp('МО спрогнозированной последовательности:');

disp(a_predict);

%Применение фильтра Калмана к оригинальной последовательности и к

%спрогнозированной

disp('ПРименение фильтра Калмана для ');

disp(filename_origin)(filename_origin,'kalmanFilter_origin',24) ;('Проверьте файл kalmanFilter_origin');('ПРименение фильтра Калмана для ');

disp(filename_predict)(filename_predict,'kalmanFilter_predict',24) ;('Проверьте файл kalmanFilter_predict');

%Вычисление ковариации от оригинальной последовательности и

%спрогнозированной=GetGaussCovariation(filename_origin,filename_predict,24) ;('Ковариация последовательностей:');

disp(cov);

7.3 Результаты

Результаты работы отобразим на рисунках 7.1 и 7.2.

Рисунок 7.1 - Первая часть результатов работы программы

Рисунок 7.2 - Вторая часть результатов работы программы

Заключение

Благодаря глубокому анализу всех основных аспектов теории вероятности и математической статистики нам удалось значительно продвинуться в понимании теоретической и прикладной составляющих самых передовых на сегодняшний день методов и концепций. Применяя эти средства как для научной, так и для промышленной разработки мы сможем достичь ещё более значительных результатов в будущем.

Литература

1.      Barker E., Kelsey J., Recommendation for Random Number Generation Using Deterministic Random Bit Generators, NIST SP800-90A, January 2012.

.        Brent R.P., "Some long-period random number generators using shifts and xors", ANZIAM Journal, 2007; 48: C188-C202.

.        Gentle J.E. (2003), Random Number Generation and Monte Carlo Methods, Springer.

4.      Hörmann W., Leydold J., Derflinger G. (2004, 2011), Automatic Nonuniform Random Variate Generation, Springer-Verlag.

.        Knuth D.E.. The Art of Computer Programming, Volume 2: Seminumerical Algorithms, Third Edition. Addison-Wesley, 1997. ISBN 0-201-89684-2. Chapter 3.

.        Luby M., Pseudorandomness and Cryptographic Applications, Princeton Univ Press, 1996.

.        Matthews R., "Maximally Periodic Reciprocals", Bulletin of the Institute of Mathematics and its Applications, 28: 147-148, 1992.

.        von Neumann J., "Various techniques used in connection with random digits," in A.S. Householder, G.E. Forsythe, and H.H. Germond, eds., Monte Carlo Method, National Bureau of Standards Applied Mathematics Series, 12 (Washington, D.C.: U.S. Government Printing Office, 1951): 36-38.

.        Peterson, Ivars (1997). The Jungles of Randomness : a mathematical safari. New York: John Wiley & Sons. ISBN 0-471-16449-6.

10.    Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. (2007), Numerical Recipes (Cambridge University Press).

.        Viega J., "Practical Random Number Generation in Software", in Proc. 19th Annual Computer Security Applications Conference, Dec. 2003.

12.    Волков И.К., Зуев С.М., Цветкова Г.М. Случайные процессы: Учеб. для вузов / Под ред. B.C. Зарубина, А.П. Крищенко. - М.: Изд-во МГТУ им. Н.Э. Баумана, 1999. - 448с.

.        ГОСТ 21878-76: Случайные процессы и динамические системы.

.        Левин Б.Р. Теоретические основы статистической радиотехники. - 3-е над.., перераб. и доп. - М.: Радио и связь, 1989. - 656с.

.        Марков А.А. Элементы математической логики. - М.: Изд-во МГУ, 1984.