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

Размещено на http://www.allbest.ru

Размещено на http://www.allbest.ru

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Пусть дано вероятностное пространство . Параметризованное семейство Размещено на http://www.allbest.ru

Размещено на http://www.allbest.ru

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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. Результаты испытаний бросания монеты

Количество опытов N	1	2	3	4	5	6	7	8	9	10	11	12	13	14
Значение счетчика выпадения орла Nо	0	0	1	1	2	3	4	…	…	…	…	…	…	…
Значение счетчика выпадения решки Nр	1	2	2	3	3	3	3	…	…	…	…	…	…	…
Частость выпадения орла Pо =Nо/N	0	0	0.33	0.25	0.4	0.5	0.57	…	…	…	…	…	…	…
Частость выпадения решки Pр =Nр/N	1	1	0.66	0.75	0.6	0.5	0.43	…	…	…	…	…	…	…

Будем подсчитывать частость выпадения орла как отношение количества случаев выпадения орла к общему числу наблюдений. Посмотрим в таблице 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 при вычислении вероятности выпадения орла

Опыт	Nкрэ
1	288
2	95
3	50
4	29
5	113
6	210
7	30
8	42
9	39
10	48
Среднее Nкр.э	94

Таким образом, проведя 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 при вычислении вероятности выпадения орла

Доверительная вероятность QF	Коэффициент Лапласа k(QF)	Требуемое число опытов Nкрт = k(QF) · p · (1 - p)/е2
0.90	2,72	68
0,95	3,84	96
0,99	6,66	167

Как видим, полученная нами оценка длины реализации, равная 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

Точность е	Критическое число экспериментов Nкрт
0,1	96
0,01	9600
0,001	960000

Построим по таблице 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. Решение задачи методом Монте-Карло

Номер точки	R	G	X	Y	Точка (X; Y) попала в прямоугольник?	Точка (X; Y) попала в пятиугольник?
1	0,8109	0,3557	8,109	7,114	Да	Да
2	0,0333	0,5370	0,333	10,740	Да	Нет
3	0,1958	0,2748	1,958	5,496	Да	Да
4	0,6982	0,1652	6,982	3,304	Да	Да
5	0,9499	0,1090	9,499	2,180	Да	Нет
6	0,7644	0,2194	7,644	4,388	Да	Да
7	0,8395	0,4510	8,395	9,020	Да	Да
8	0,0415	0,6855	0,415	13,710	Да	Нет
9	0,5997	0,1140	5,997	2,280	Да	Да
10	0,9595	0,9595	9,595	19,190	Да	Нет
Всего:	10	6

Статистическая гипотеза заключается в том, что количество точек, попавших в контур фигуры, пропорционально площади фигуры: 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. Изучение алгоритмов компьютерного моделирования случайных последовательностей

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 - Диаграмма получения случайных чисел аппаратным методом

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

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

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

	Источник энтропии	ГПСЧ	Достоинства	Недостатки
/dev/random в UNIX/Linux	Счётчик тактов процессора, однако собирается только во время аппаратных прерываний	LFSR, с хешированием выхода через SHA-1	Есть во всех Unix, надёжный источник энтропии	Очень долго «нагревается», может надолго «застревать», либо работает как ГПСЧ (/dev/urandom)
Yarrow от Брюса Шнайера	Традиционные методы	AES-256 и SHA-1 маленького внутреннего состояния	Гибкий криптостойкий дизайн	Медленный
Microsoft Crypto API	Текущее время, размер жёсткого диска, размер свободной памяти, номер процесса и NETBIOS-имя компьютера	MD5-хеш внутреннего состояния размером в 128 бит	Встроен в Windows, не «застревает»	Сильно зависит от используемого криптопровайдера (CSP).
Java Secure Random	Взаимодействие между потоками	SHA-1-хеш внутреннего состояния (1024 бит)	Большое внутреннее состояние	Медленный сбор энтропии
Chaos от Ruptor	Счётчик тактов процессора, собирается непрерывно	Хеширование 4096-битового внутреннего состояния на основе нелинейного варианта Marsaglia-генератора	Пока самый быстрый из всех, большое внутреннее состояние, не «застревает»	Оригинальная разработка, свойства приведены только по утверждению автора
RRAND от Ruptor	Счётчик тактов процессора	Зашифровывание внутреннего состояния поточным шифром EnRUPT в authenticated encryption режиме (aeRUPT)	Очень быстр, внутреннее состояние произвольного размера по выбору, не «застревает»	Оригинальная разработка, свойства приведены только по утверждению автора. Шифр EnRUPT не является криптостойким
RdRand от Intel	Шумы токов	Построение ПСЧ на основе «случайного» битового считывания значений от токов	Очень быстр, не «застревает»	Оригинальная разработка, свойства приведены только по утверждению статьи из habrahabr
ГПСЧ Stratosphera от ORION	Счетчик тактов процессора, собирается непрерывно (также используется соль в виде случайно выбранного целого	Построение ПСЧ на основе алгоритма от Intel с многоразовой инициализацией и сдвигом	Достаточно быстр, не «застревает», проходит все тесты DIEHARD	Оригинальная разработка, свойства приведены только исходя из информации на сайте oriondevteam.com.

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

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

Случайные цифры	Равномерно распределенные от 0 до 1 случайные числа
9	2	9	2	0	4	2	6	0,929
9	5	7	3	4	9	0	3	0,204
5	9	1	6	6	5	7	6	0,269

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

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

ri + 1 = f(ri).

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

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

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

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

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

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

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

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

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

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

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

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

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

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.6).

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

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

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

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

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

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

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

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. Область распределения может быть ограничена любым интервалом.

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

Если же системно (применяя достижения алгебры и численных методов, опирающиеся на твёрдую поверхность кибернетики) подходить к каждому конкретному распределению, то на сегодняшний день можно отметить работы Александра Самарина (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), и выполняете следующий алгоритм:

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

2. Если 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. Создание программ моделирования случайных последовательностей

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Содержимое файла 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); % марковская последовательность

end

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

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

disp(answer);

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

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



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

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

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

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

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

г = n(t)

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

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

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

t2 > 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'), то двумерная плотность вероятности равна:

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

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

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

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

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