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

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

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

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

function GaussScalarSequense(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) ;

disp(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(), которая возвращает сумму чисел. Сделаем её нашей отправной точкой, и придём к успеху (путём деления суммы чисел на их количество (подходящий тип оценки матожидания)).

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

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

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

Содержимое файла 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

summ=0;

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

summ=sum(values);

fclose(fid);

a=summ/N;

end

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

GaussScalarSequense('Gauss1',0,1,1000) ;

answer=GetMiddleOfGauss('Gauss1',1000) ;

disp('Middle of gausse sequense:');

disp(answer);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

function covar= 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;

end

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

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

GaussScalarSequense('Gauss2',0,2,100) ;

answer=GetGaussCovariation('Gauss1','Gauss2',100) ;

disp('Covariation of gausse sequenseS:');

disp(answer);

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

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



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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5.2 Проектирование

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

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

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

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

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

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

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

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

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

5.3 Разработка

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

clear; clc;

% <Date> - <Hour> - <Value>

load prices PRICES_EUR; % Загрузка файла prices.mat

TimeSeries = PRICES_EUR;

% Задача прогнозирования

T = datenum('01.09.2012 23:00:00', 'dd.mm.yyyy hh:MM:ss'); % Момент прогноза

P = 24; % Горизонт прогнозирования (Forecast)

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

M = 144;

Step = 24;

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

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

if 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;

while 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);

MSP = Likeness(IndexLikeness(1),1);

fprintf(['Максимальное подобие = ', 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);

E = X(:,2);

Xn = X'* X;

Yn = X'* Y;

invX = inv(Xn);

A = 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). Используется для уменьшения размерности задачи;

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

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

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

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

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

6.2 Проектирование

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

У нас в распоряжении есть встроенные возможности развитого САПР 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('Гауссова последователность, подвергнутая фильтру Калмана:');

disp(answer);

%Параметры

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

6.4 Результаты

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

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

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

7.1 Проектирование

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

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

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

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

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

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

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

7.2 Разработка

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

clear; clc;

% <Date> - <Hour> - <Value>

load prices PRICES_EUR; % Загрузка файла prices.mat

TimeSeries = PRICES_EUR;

% Задача прогнозирования

T = datenum('01.09.2012 23:00:00', 'dd.mm.yyyy hh:MM:ss'); % Момент прогноза

P = 24; % Горизонт прогнозирования (Forecast)

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

M = 144;

Step = 24;

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

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

if 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;

while 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);

MSP = 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);

E = X(:,2);

Xn = X'* X;

Yn = X'* Y;

invX = inv(Xn);

A = invX * Yn; % Искомые коэффициенты alpha1 и alpha0 являются значениями матрица A

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

X = HistBaseData(:,3);

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

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

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

disp(Forecast)

% 1) Определяем фактические значения

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

Fact = TimeSeries([Index : Index+P-1],3);

disp('Реальная последовательность:')

disp(Fact)

filename_origin='original_sequence';

filename_predict='predict_sequence';

fid_origin = fopen(filename_origin, 'w');

fid_predict = fopen(filename_predict, 'w');

if fid_origin == -1

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

end

if fid_predict == -1

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

end

for i=1:24

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

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

end

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

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

a_origin=GetMiddleOfGauss(filename_origin,24) ;

a_predict=GetMiddleOfGauss(filename_predict,24) ;

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

disp(a_origin);

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

disp(a_predict);

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

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

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

disp(filename_origin)

KalmanFilter(filename_origin,'kalmanFilter_origin',24) ;

disp('Проверьте файл kalmanFilter_origin');

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

disp(filename_predict)

KalmanFilter(filename_predict,'kalmanFilter_predict',24) ;

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

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

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

cov=GetGaussCovariation(filename_origin,filename_predict,24) ;

disp('Ковариация последовательностей:');

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.

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

3. 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.

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

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

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

8. 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.

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

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

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

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

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

15. Марков А.А., Нагорный Н.М. Теория алгорифмов, изд. 2. - М.: ФАЗИС, 1996.

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

Размещено на Allbest.ru