Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- funcprot(0);
- clc;
- FILE_NAME = 'C:\Users\MStul\Desktop\Курсовая работа\Array_1.txt';
- FILE_CORR = 'C:\Users\MStul\Desktop\Курсовая работа\R(10).txt';
- FLOAT_FORMAT = '%16.4f'; //формат вывода дробных чисел
- INT_FORMAT = '%d'; //формат вывода целых чисел
- //ввод данных и задание глобальных параметров, которые используются во всей программе без изменения
- x = fscanfMat(FILE_NAME); //вектор исходных данных
- n = length(x); //размер вектора
- mX = mean(x); //выборочное среднее
- //оценка моментных функций
- //дисперсия
- function D=dispersion(samp)
- D=0.0;
- mS=mean(samp);
- k=size(samp,1);
- for i=1:k
- D=D+(samp(i)-mS)*(samp(i)-mS);
- end;
- D=D/(k-1);
- endfunction
- //корреляционная функция R(k)
- function R=correlation(samp, k)
- if (k<0) then
- k=k*(-1);
- end
- R=0.0;
- r=size(samp, 1);
- mS=mean(samp);
- for i=1:(r-k)
- R=R+(samp(i)-mS)*(samp(i+k)-mS);
- end;
- R=R/(r-k-1);
- endfunction
- //вычисление нормированной корреляционной функции
- function nR=normCorrelation(samp, k)
- nR=correlation(samp,k)/dispersion(samp);
- endfunction
- //радиус корреляции
- function T=correlationRadius()
- coefficient=0.01;
- T = coefficient * length(x) - 2;
- em1 = exp(-1);
- while (T >= 0) & (abs(normCorrelation(x,T)) < em1),
- T = T - 1;
- end;
- T = T + 1;
- endfunction
- //построение графика случайного процесса
- function processPlot()
- t=0:120;
- y = zeros(length(t), 1);
- dX = zeros(length(t), 1);
- dX_1 = zeros(length(t), 1);
- mXPlot = zeros(length(t), 1);
- F=figure('figure_name','Случайный процесс', 'background',8);
- for i=1:length(t)
- y(i)=x(i);
- dX(i)=mX-dispersion(x)^(0.5);
- dX_1(i)=mX+dispersion(x)^(0.5);
- mXPlot(i)=mX;
- end;
- plot2d(t,y, axesflag = 1, style = 2);
- plot2d(t,mXPlot, axesflag = 1, style = 5);
- plot2d(t,dX, axesflag = 1, style = 3);
- plot2d(t,dX_1, axesflag = 1, style = 3);
- xtitle("","Индекс","Случайный процесс");
- xgrid();
- legend('Исходный процесс', 'Выборочное среднее', 'СКО',5);
- endfunction
- //построение графика нормированной корреляционной функции
- function normCorrPlot()
- t = [0:10];
- y = zeros(length(t), 1);
- exp_p = zeros(length(t), 1);
- exp_m = zeros(length(t), 1);
- for i = 1:length(t),
- y(i) = normCorrelation(x,t(i));
- exp_p(i) = exp(-1);
- exp_m(i) = -exp(-1);
- end;
- d = [0:1];
- interval = [-0.3 -0.3];
- F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
- plot2d(t,y, axesflag = 1,style = 2);
- plot2d(t,exp_p, axesflag = 1, style = 3);
- plot2d(t,exp_m, axesflag = 1, style = 3);
- plot2d(d, interval, alexflag = 1, style = 5);
- xtitle("","Индекс","Нормированная корреляционная функция");
- xgrid();
- legend('НКФ', '+1/e -1/e', 'Интервал корреляции');
- endfunction;
- //вывод на экран всех расчётов
- function Program1()
- printf("ЗАДАНИЕ 1------------------------------\n");
- processPlot();
- printf("Выборочное среднее: " + FLOAT_FORMAT + "\n", mX);
- dX=dispersion(x);
- printf("Выборочная дисперсия: " + FLOAT_FORMAT + "\n", dX);
- printf("CКО:" + FLOAT_FORMAT + "\n", sqrt(dX));
- for k=0:10
- printf("Корреляционная функция для k=" + INT_FORMAT + FLOAT_FORMAT + "\n", k, correlation(x,k));
- printf("Нормированная корреляционная функция для k=" + INT_FORMAT + FLOAT_FORMAT + "\n", k, normCorrelation(x,k));
- end;
- printf("Радиус корреляции: " + FLOAT_FORMAT + "\n", correlationRadius());
- normCorrPlot();
- endfunction
- //Вспомогательные функции для выполнения следующих заданий
- //получение вектора корреляционных функций от R(0) до R(k)
- function RVector=correlationVector(samp, k)
- RVector=zeros(1,k+1);
- for i =0:k
- RVector(i+1)=correlation(samp,i)
- end
- endfunction
- //получение вектора нормированных корреляционных функций от R(0) до R(k)
- function nRVector=normCorrelationVector(samp, k)
- nRVector=zeros(1,k+1);
- for i =0:k
- nRVector(i+1)=normCorrelation(samp,i)
- end
- endfunction
- //2. Построение моделей АР(M) и выбор наилучшей модели
- //Нахождение коэффициентов betta и alpha
- function betas=foundBetas(M, N)
- betas = zeros(1, M);
- A = zeros(M, M);
- b = zeros(1, M);
- for i = 1 : M
- for j = 1 : M
- A(i, j) = correlation(x,N + i - j);
- end
- b(i) = -correlation(x,N + i);
- end
- betas = linsolve(A, b');
- endfunction
- //Нахождение коэффициента alpha
- function alpha=calcAlpha(betas)
- alpha = correlation(x,0);
- M=size(betas,1);
- for i = 1 : M
- alpha = alpha - betas(i) * correlation(x,i);
- end
- alpha = sqrt(alpha);
- endfunction
- //проверка построенной модели на устойчивость
- function result=checkStability(betas)
- sizeB = length(betas);
- if (sizeB == 0) then result = %T;
- elseif ((sizeB == 1) & (abs(betas(1)) < 1)) then result = %T;
- elseif ((sizeB == 2) & (abs(betas(2)) < 1) & (abs(betas(1)) < 1 - betas(2))) then result = %T;
- elseif ((sizeB == 3) & (abs(betas(3)) < 1)) & (abs(betas(1) + betas(3)) < 1 - betas(2)) & (abs(betas(2) + betas(1) * betas(3)) < abs(1 - betas(3)^2)) then result = %T;
- else
- result = %F;
- end
- endfunction
- //вычисление теоретической корреляционной функции R(k)
- function tR=theorCorrelation(RX, alphas, betas, k)
- nm = length(alphas)+length(betas)-1;
- k = abs(k);
- if (k > nm) then
- tR = 0;
- M = length(betas);
- for j = 1 : M,
- tR = tR + betas(j) * theorCorrelation(RX,alphas,betas, k - j);
- end;
- else
- tR = RX(k+1);
- end
- endfunction;
- //Вычисление нормальной теоретической корреляционной функции r(k)
- //(используется так же для пунктов 3,4,5,6)
- function tnR=theorNormCorrelation(RX, alphas, betas, k)
- tnR=zeros(k+1);
- for i = 0:k
- tnR(i+1) = theorCorrelation(RX,alphas,betas,i);
- end;
- R=theorCorrelation(RX,alphas,betas,0);
- tnR = tnR./R;
- endfunction
- //вычисление среднего квадратичного отклонения
- //(используется так же для пунктов 3,4,5,6)
- function epsilon=epsilonForNCF(r, tr)
- epsilon = 0;
- m = min(length(r), length(tr));
- for j = 1 : m,
- epsilon = epsilon + (r(j) - tr(j))^2;
- end;
- endfunction;
- //Вывод на экран всех расчётов
- function minM=Program2(M)
- printf("ЗАДАНИЕ 2------------------------------\n");
- qE=zeros(M);
- for i =1:M
- betas = foundBetas(i, 0);
- alpha=calcAlpha(betas)
- printf("Коэффициенты betta для M=" + INT_FORMAT + ":\n", i);
- if (checkStability(betas)) then printf(FLOAT_FORMAT + "\n", betas);
- else printf("Система нестабильна");
- end
- printf("Коэффициент alpha для M=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, alpha);
- t=normCorrelationVector(x,10);
- tr=theorNormCorrelation(correlationVector(x,10),alpha, betas, 10);
- qE(i)=epsilonForNCF(t,tr);
- printf("Среднее квадратичное отклонение для M=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, qE(i));
- end
- minM=qE(1);
- for i=2:M
- if minM>qE(i) then
- minM=qE(i);
- minM=i;
- end
- end
- printf("Лучшей моделью является модель АР(" + INT_FORMAT + ")\n", minM);
- endfunction
- //3. Построение моделей СС(N) и выбор наилучшей модели
- //создание нелинейной системы для нахождения коэффициентов alpha
- function y=createNoLinearSystem(alphas, N, R)
- y = zeros(N + 1, 1);
- for i = 1 : N + 1
- for j = 1 : (N - i + 2)
- y(i) = y(i) + alphas(j) * alphas(j + i - 1);
- end
- y(i) = y(i) - R(i);
- end
- endfunction
- //Нахождение коэффициентов alpha
- function alphas=calcAlphas(N)
- R=correlationVector(x,10);
- alphas = zeros(N + 1, 1);
- alphas = fsolve(alphas, list(createNoLinearSystem, N, R));
- endfunction
- //Вывод на экран всех расчётов
- function minN=Program3(N)
- printf("ЗАДАНИЕ 3------------------------------\n");
- qE=zeros(N+1);
- for i =0:N
- alphas = calcAlphas(i);
- printf("Коэффициенты alpha для N=" + INT_FORMAT + ":\n", i);
- printf(FLOAT_FORMAT + "\n", alphas);
- t=normCorrelationVector(x,10);
- tr=theorNormCorrelation(correlationVector(x,10),alphas, [], 10);
- qE(i+1)=epsilonForNCF(t,tr);
- printf("Среднее квадратичное отклонение для N=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, qE(i+1));
- end
- minN=qE(1);
- for i=2:N
- if minN>qE(i) then
- minN=qE(i);
- minN=i;
- end
- end
- printf("Лучшей моделью является модель CC(" + INT_FORMAT + ")\n", minN);
- endfunction
- //4 задание
- R = fscanfMat(FILE_CORR); //вектор корреляционной функции R(10)
- function f = arma11(p)
- a0 = p(1); a1 = p(2); b1 = p(3);
- x0 = p(4); x1 = p(5);
- f(1) = b1*R(2) + a0*x0 + a1*x1 - R(1);
- f(2) = b1*R(1) + a1*x0 - R(2);
- f(3) = b1*R(2) - R(3);
- f(4) = a0 - x0;
- f(5) = b1*x0 + a1 - x1;
- end;
- function f = arma12(p)
- a0 = p(1); a1 = p(2); a2 = p(3); b1 = p(4); x0 = p(5); x1 = p(6); x2 = p(7);
- f(1) = b1*R(2) + a0*x0 + a1*x1 + a2*x2 - R(1);
- f(2) = b1*R(1) + a1*x0 + a2*x1 - R(2);
- f(3) = b1*R(2) + a2*x0 - R(3);
- f(4) = b1*R(3) - R(4);
- f(5) = a0 - x0;
- f(6) = b1*x0 + a1 - x1;
- f(7) = b1*x1 + a2 - x2;
- end;
- function f = arma13(p)
- a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4);
- b1 = p(5); x0 = p(6); x1 = p(7); x2 = p(8); x3 = p(9);
- f(1) = b1*R(2) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
- f(2) = b1*R(1) + a1*x0 + a2*x1 + a3*x2 - R(2);
- f(3) = b1*R(2) + a2*x0 + a3*x1 - R(3);
- f(4) = b1*R(3) + a3*x0 - R(4);
- f(5) = b1*R(4) - R(5);
- f(6) = a0 - x0;
- f(7) = b1*x0 + a1 - x1;
- f(8) = b1*x1 + a2 - x2;
- f(9) = b1*x1 + a3 - x3;
- end;
- function f = arma21(p)
- a0 = p(1); a1 = p(2); b1 = p(3); b2 = p(4);
- x0 = p(5); x1 = p(6);
- f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 - R(1);
- f(2) = b1*R(1) + b2*R(2) + a1*x0 - R(2);
- f(3) = b1*R(2) + b2*R(1) - R(3);
- f(4) = b1*R(3) + b2*R(2) - R(4);
- f(5) = a0 - x0;
- f(6) = b1*x0 + a1 - x1;
- end;
- function f = arma22(p)
- a0 = p(1); a1 = p(2); a2 = p(3); b1 = p(4);
- b2 = p(5); x0 = p(6); x1 = p(7); x2 = p(8);
- f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 + a2*x2 - R(1);
- f(2) = b1*R(1) + b2*R(2) + a1*x0 + a2*x1 - R(2);
- f(3) = b1*R(2) + b2*R(1) + a2*x0 - R(3);
- f(4) = b1*R(3) + b2*R(2) - R(4);
- f(5) = b1*R(4) + b2*R(3) - R(5);
- f(6) = a0 - x0;
- f(7) = b1*x0 + a1 - x1;
- f(8) = b1*x1 + b2*x0 + a2 - x2;
- end;
- function f = arma23(p)
- a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4); b1 = p(5);
- b2 = p(6); x0 = p(7); x1 = p(8); x2 = p(9); x3 = p(10);
- f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
- f(2) = b1*R(1) + b2*R(2) + a1*x0 + a2*x1 + a3*x2 - R(2);
- f(3) = b1*R(2) + b2*R(1) + a2*x0 + a3*x1 - R(3);
- f(4) = b1*R(3) + b2*R(2) + a3*x0 - R(4);
- f(5) = b1*R(4) + b2*R(3) - R(5);
- f(6) = b1*R(5) + b2*R(4) - R(6);
- f(7) = a0 - x0;
- f(8) = b1*x0 + a1 - x1;
- f(9) = b1*x1 + b2*x0 + a2 - x2;
- f(10) = b1*x2 + b2*x1 + a3 - x3;
- end;
- function f = arma31(p)
- a0 = p(1); a1 = p(2); b1 = p(3); b2 = p(4);
- b3 = p(5); x0 = p(6); x1 = p(7);
- f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 - R(1);
- f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 - R(2);
- f(3) = b1*R(2) + b2*R(1) + b3*R(2) - R(3);
- f(4) = b1*R(3) + b2*R(2) + b3*R(1) - R(4);
- f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
- f(6) = a0 - x0;
- f(7) = b1*x0 + a1 - x1;
- end;
- function f = arma32(p)
- a0 = p(1); a1 = p(2); a2 = p(3); b1 = p(4); b2 = p(5);
- b3 = p(6); x0 = p(7); x1 = p(8); x2 = p(9);
- f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 + a2*x2 - R(1);
- f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 + a2*x1 - R(2);
- f(3) = b1*R(2) + b2*R(1) + b3*R(2) + a2*x0 - R(3);
- f(4) = b1*R(3) + b2*R(2) + b3*R(1) - R(4);
- f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
- f(6) = b1*R(5) + b2*R(4) + b3*R(3) - R(6);
- f(7) = a0 - x0;
- f(8) = b1*x0 + a1 - x1;
- f(9) = b1*x1 + b2*x0 + a2 - x2;
- end;
- function f = arma33(p)
- a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4); b1 = p(5); b2 = p(6);
- b3 = p(7); x0 = p(8); x1 = p(9); x2 = p(10); x3 = p(11);
- f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
- f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 + a2*x1 + a3*x2 - R(2);
- f(3) = b1*R(2) + b2*R(1) + b3*R(2) + a2*x0 + a3*x1 - R(3);
- f(4) = b1*R(3) + b2*R(2) + b3*R(1) + a3*x0 - R(4);
- f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
- f(6) = b1*R(5) + b2*R(4) + b3*R(3) - R(6);
- f(7) = b1*R(6) + b2*R(5) + b3*R(4) - R(7);
- f(8) = a0 - x0;
- f(9) = b1*x0 + a1 - x1;
- f(10) = b1*x1 + b2*x0 + a2 - x2;
- f(11) = b1*x2 + b2*x1 + b3*x0 + a3 - x3;
- end;
- function [koef, bestM, bestN] = arma_coefficients(R)
- r = R / R(1);
- qEr = zeros(1, 9);
- printf('\nCreate and research ARMA model:\n');
- printf('%8s %8s %8s %8s %8s %8s %8s ', 'a0', 'a1', 'a2', 'a3', 'b1', 'b2', 'b3');
- printf('\n');
- [answer,fun,info(1)] = fsolve([sqrt(R(1));0;0;0;0], arma11);
- if (info(1) <> 4 && norm(fun) < 1e-4)
- koef(1, :) = [answer(1), answer(2), 0, 0, answer(3), 0, 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(1, :));
- alphas = koef(1, 1:2)
- betas = koef(1, 5);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(1) = epsilonForNCF(r, tr);
- else
- info(1) = 0;
- printf('ARMA(1,1) does not exists');
- koef(1, :) = zeros(1,7);
- qEr(1) = %inf;
- end
- printf('\n');
- [answer,fun,info(2)] = fsolve([sqrt(R(1));0;0;0;0;0;0], arma12);
- if (info(2) <> 4 && norm(fun) < 1e-4)
- koef(2, :) = [answer(1), answer(2), answer(3), 0, answer(4), 0, 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(2, :));
- alphas = koef(2, 1:3)
- betas = koef(2, 5);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(2) = epsilonForNCF(r, tr);
- else
- info(2) = 0;
- printf('ARMA(1,2) does not exists');
- koef(2, :) = zeros(1,7);
- qEr(2) = %inf;
- end
- printf('\n');
- [answer,fun,info(3)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0], arma13);
- if (info(3) <> 4 && norm(fun) < 1e-4)
- koef(3, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), 0, 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(3, :));
- alphas = koef(3, 1:4)
- betas = koef(3, 5);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(3) = epsilonForNCF(r, tr);
- else
- info(3) = 0;
- printf('ARMA(1,3) does not exists');
- koef(3, :) = zeros(1,7);
- qEr(3) = %inf;
- end
- printf('\n');
- [answer,fun,info(4)] = fsolve([sqrt(R(1));0;0;0;0;0], arma21);
- if (info(4) <> 4 && norm(fun) < 1e-4)
- koef(4, :) = [answer(1), answer(2), 0, 0, answer(3), answer(4), 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(4, :));
- alphas = koef(4, 1:2);
- betas = koef(4, 5:6);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(4) = epsilonForNCF(r, tr);
- else
- info(4) = 0;
- printf('ARMA(2,1) does not exists');
- koef(4, :) = zeros(1,7);
- qEr(4) = %inf;
- end
- printf('\n');
- [answer,fun,info(5)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0], arma22);
- if (info(5) <> 4 && norm(fun) < 1e-4)
- koef(5, :) = [answer(1), answer(2), answer(3), 0, answer(4), answer(5), 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(5, :));
- alphas = koef(5, 1:3);
- betas = koef(5, 5:6);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(5) = epsilonForNCF(r, tr);
- else
- info(5) = 0;
- printf('ARMA(2,2) does not exists');
- koef(5, :) = zeros(1,7);
- qEr(5) = %inf;
- end
- printf('\n');
- [answer,fun,info(6)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0;0], arma23);
- if (info(6) <> 4 && norm(fun) < 1e-4)
- koef(6, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), answer(6), 0];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(6, :));
- alphas = koef(6, 1:4);
- betas = koef(6, 5:6);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(6) = epsilonForNCF(r, tr);
- else
- info(6) = 0;
- printf('ARMA(2,3) does not exists');
- koef(6, :) = zeros(1,7);
- qEr(6) = %inf;
- end
- printf('\n');
- [answer,fun,info(7)] = fsolve([sqrt(R(1));0;0;0;0;0;0], arma31);
- if (info(7) <> 4 && norm(fun) < 1e-4)
- koef(7, :) = [answer(1), answer(2), 0, 0, answer(3), answer(4), answer(5)];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(7, :));
- alphas = koef(7, 1:2);
- betas = koef(7, 5:7);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(7) = epsilonForNCF(r, tr);
- else
- info(7) = 0;
- printf('ARMA(3,1) does not exists');
- koef(7, :) = zeros(1,7);
- qEr(7) = %inf;
- end
- printf('\n');
- [answer,fun,info(8)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0], arma32);
- if (info(8) <> 4 && norm(fun) < 1e-4)
- koef(8, :) = [answer(1), answer(2), answer(3), 0, answer(4), answer(5), answer(6)];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(8, :));
- alphas = koef(8, 1:3);
- betas = koef(8, 5:7);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- qEr(8) = epsilonForNCF(r, tr);
- else
- info(8) = 0;
- printf('ARMA(3,2) does not exists');
- koef(8, :) = zeros(1,7);
- qEr(8) = %inf;
- end
- printf('\n');
- [answer,fun,info(9)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0;0;0], arma33);
- if (info(9) <> 4 && norm(fun) < 1e-4)
- koef(9, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), answer(6), answer(7)];
- printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(9, :));
- alphas = koef(9, 1:4);
- betas = koef(9, 5:7);
- tr = theorNormCorrelation(R, alphas, betas, 10);
- //disp(r);
- //disp(tr);
- qEr(9) = epsilonForNCF(r, tr);
- else
- info(9) = 0;
- printf('ARMA(3,3) does not exists');
- koef(9, :) = zeros(1,7);
- qEr(9) = %inf;
- end
- printf('\n');
- for i=1:9
- if (qEr(i) == %inf)
- printf('Inf\n');
- else
- printf('%.4f\n', qEr(i));
- end
- end
- printf('\n');
- //выберем модель с минимальным значением qEr - это будет наилучшая модель
- minqEr = min(qEr);
- best = -1;
- for i = 1:9
- if (qEr(i) == minqEr)
- best = i;
- end
- end
- bestM = floor((best - 1)/3) + 1; //floor - функция округления вниз
- bestN = modulo(best - 1, 3) + 1;
- koef = koef(best,:);
- endfunction
- [coefARMA, armaM, armaN] = arma_coefficients(R);
- function res = createProcessValues(coefProcess, armaM, armaN)
- r=rand(1, 11000,"normal"); //генерируем 11000 значений, распределенных по нормальному закону
- res = zeros(1, 10000);
- ksi = r(1001:11000);
- for i = 1:10000
- j = 1;
- //этта
- while (j <= armaM && i > j)
- res(i) = res(i) + coefARMA(4+j) * res(i-j); //coefARMA(4+j); берем 4, т.к. беты записаны 5, 6, 7 элементом
- j = j + 1;
- end
- d = 1;
- //кси
- while (d <= armaN && i >= d)
- res(i) = res(i) + coefARMA(d) * ksi(i-d+1); //coefARMA(d); берем d сразу, т.к альфы записаны под 1, 2, 3, 4
- d = d + 1;
- end
- //математическое ожидание
- sumTemp = 0;
- for k = 1:armaM
- sumTemp = sumTemp + coefARMA(4+k);
- end
- res(i) = res(i) + mX * (1 - sumTemp);
- end
- endfunction
- resARMA = createProcessValues(coefARMA, armaM, armaN);
- //запишем коэффициенты наилучшей модели в переменную
- minAR = Program2(3);
- minMA = Program3(3);
- betasAR = foundBetas(minAR, 0);
- alphaAR=calcAlpha(betasAR);
- coefAR = [alphaAR(1), 0, 0, 0, betasAR(1), betasAR(2), betasAR(3)];
- alphaMA = calcAlphas(minMA);//для СС
- coefMA = [alphaMA(1), alphaMA(2), alphaMA(3), alphaMA(4), 0, 0, 0];
- resAR = createProcessValues(coefAR, minAR, 0);
- resMA = createProcessValues(coefMA, 0, minMA);
- //построение графиков для 5 задания
- //график для АРСС
- function processPlotForARMA()
- t = [0:10];
- alpha = coefARMA(1:4);
- betas = coefARMA(5:7);
- normCorr_x=normCorrelationVector(x,10);
- theorNormCorr=theorNormCorrelation(correlationVector(x,10),alpha, betas, 10);
- normCorr_res=normCorrelationVector(resARMA',10);
- F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
- plot2d(t,normCorr_x, axesflag = 1,style = 5);
- plot2d(t,theorNormCorr, axesflag = 1,style = 3);
- plot2d(t,normCorr_res, axesflag = 1,style = 2);
- xtitle("","Индекс","Нормированная корреляционная функция");
- xgrid();
- legend('Процесс', 'Теоретическая АРСС(3,3)', 'Смоделированный процесс');
- endfunction;
- processPlotForARMA();
- //график для АР
- function processPlotForAR()
- t = [0:10];
- normCorr_x=normCorrelationVector(x,10);
- theorNormCorr=theorNormCorrelation(correlationVector(x,10),alphaAR, betasAR, 10);
- normCorr_res=normCorrelationVector(resAR',10);
- F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
- plot2d(t,normCorr_x, axesflag = 1,style = 5);
- plot2d(t,theorNormCorr, axesflag = 1,style = 3);
- plot2d(t,normCorr_res, axesflag = 1,style = 2);
- xtitle("","Индекс","Нормированная корреляционная функция");
- xgrid();
- legend('Процесс', 'Теоретическая АР(3)', 'Смоделированный процесс');
- endfunction;
- processPlotForAR();
- //график для СС
- function processPlotForMA()
- t = [0:10];
- normCorr_x=normCorrelationVector(x,10);
- theorNormCorr=theorNormCorrelation(correlationVector(x,10),alphaMA, [], 10);
- normCorr_res=normCorrelationVector(resMA',10);
- F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
- plot2d(t,normCorr_x, axesflag = 1,style = 5);
- plot2d(t,theorNormCorr, axesflag = 1,style = 3);
- plot2d(t,normCorr_res, axesflag = 1,style = 2);
- xtitle("","Индекс","Нормированная корреляционная функция");
- xgrid();
- legend('Процесс', 'Теоретическая СС(3)', 'Смоделированный процесс');
- endfunction;
- processPlotForMA();
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement