Advertisement
Guest User

kursach_1

a guest
Apr 20th, 2019
171
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Scilab 23.52 KB | None | 0 0
  1. funcprot(0);
  2. clc;
  3. FILE_NAME = 'C:\Users\MStul\Desktop\Курсовая работа\Array_1.txt';
  4. FILE_CORR = 'C:\Users\MStul\Desktop\Курсовая работа\R(10).txt';
  5. FLOAT_FORMAT = '%16.4f'; //формат вывода дробных чисел
  6. INT_FORMAT = '%d';  //формат вывода целых чисел
  7.  
  8. //ввод данных и задание глобальных параметров, которые используются во всей программе без изменения
  9. x = fscanfMat(FILE_NAME); //вектор исходных данных
  10. n = length(x); //размер вектора
  11. mX = mean(x); //выборочное среднее
  12.  
  13. //оценка моментных функций
  14. //дисперсия
  15. function D=dispersion(samp)
  16.     D=0.0;
  17.     mS=mean(samp);
  18.     k=size(samp,1);
  19.     for i=1:k
  20.         D=D+(samp(i)-mS)*(samp(i)-mS);
  21.     end;
  22.     D=D/(k-1);
  23. endfunction
  24.  
  25. //корреляционная функция R(k)
  26. function R=correlation(samp, k)
  27.         if (k<0) then
  28.         k=k*(-1);
  29.     end
  30.     R=0.0;
  31.     r=size(samp, 1);
  32.     mS=mean(samp);
  33.     for i=1:(r-k)
  34.         R=R+(samp(i)-mS)*(samp(i+k)-mS);
  35.     end;
  36.     R=R/(r-k-1);
  37. endfunction
  38.  
  39. //вычисление нормированной корреляционной функции
  40. function nR=normCorrelation(samp, k)
  41.     nR=correlation(samp,k)/dispersion(samp);
  42. endfunction
  43.  
  44. //радиус корреляции
  45. function T=correlationRadius()
  46.     coefficient=0.01;
  47.     T = coefficient * length(x) - 2;
  48.     em1 = exp(-1);
  49.     while (T >= 0) & (abs(normCorrelation(x,T)) < em1),
  50.         T = T - 1;
  51.     end;
  52.     T = T + 1;
  53. endfunction
  54.  
  55. //построение графика случайного процесса
  56. function processPlot()
  57.     t=0:120;
  58.     y = zeros(length(t), 1);
  59.     dX = zeros(length(t), 1);
  60.     dX_1 = zeros(length(t), 1);
  61.     mXPlot = zeros(length(t), 1);
  62.     F=figure('figure_name','Случайный процесс', 'background',8);
  63.     for i=1:length(t)
  64.         y(i)=x(i);
  65.         dX(i)=mX-dispersion(x)^(0.5);
  66.         dX_1(i)=mX+dispersion(x)^(0.5);
  67.         mXPlot(i)=mX;
  68.     end;       
  69.     plot2d(t,y, axesflag = 1, style = 2);
  70.     plot2d(t,mXPlot, axesflag = 1, style = 5);
  71.     plot2d(t,dX, axesflag = 1, style = 3);
  72.     plot2d(t,dX_1, axesflag = 1, style = 3);
  73.     xtitle("","Индекс","Случайный процесс");
  74.     xgrid();
  75.     legend('Исходный процесс', 'Выборочное среднее', 'СКО',5);
  76. endfunction
  77.  
  78. //построение графика нормированной корреляционной функции
  79. function normCorrPlot()
  80.     t = [0:10];
  81.     y = zeros(length(t), 1);
  82.     exp_p = zeros(length(t), 1);
  83.     exp_m = zeros(length(t), 1);
  84.     for i = 1:length(t),
  85.         y(i) = normCorrelation(x,t(i));
  86.     exp_p(i) = exp(-1);
  87.     exp_m(i) = -exp(-1);
  88.     end;
  89.  
  90.     d = [0:1];
  91.     interval = [-0.3 -0.3];
  92.     F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
  93.     plot2d(t,y, axesflag = 1,style = 2);
  94.     plot2d(t,exp_p, axesflag = 1, style = 3);
  95.     plot2d(t,exp_m, axesflag = 1, style = 3);    
  96.     plot2d(d, interval, alexflag = 1, style = 5);
  97.     xtitle("","Индекс","Нормированная корреляционная функция");
  98.     xgrid();
  99.     legend('НКФ', '+1/e -1/e', 'Интервал корреляции');
  100. endfunction;
  101.  
  102. //вывод на экран всех расчётов
  103. function Program1()
  104.     printf("ЗАДАНИЕ 1------------------------------\n");
  105.     processPlot();
  106.     printf("Выборочное среднее: " + FLOAT_FORMAT + "\n", mX);
  107.     dX=dispersion(x);
  108.     printf("Выборочная дисперсия: " + FLOAT_FORMAT + "\n", dX);
  109.     printf("CКО:" + FLOAT_FORMAT + "\n", sqrt(dX));
  110.     for k=0:10
  111.         printf("Корреляционная функция для k=" + INT_FORMAT + FLOAT_FORMAT + "\n", k, correlation(x,k));
  112.         printf("Нормированная корреляционная функция для k=" + INT_FORMAT + FLOAT_FORMAT + "\n", k, normCorrelation(x,k));
  113.     end;
  114.     printf("Радиус корреляции: " + FLOAT_FORMAT + "\n", correlationRadius());
  115.     normCorrPlot();
  116. endfunction
  117.  
  118. //Вспомогательные функции для выполнения следующих заданий
  119. //получение вектора корреляционных функций от R(0) до R(k)
  120. function RVector=correlationVector(samp, k)
  121.     RVector=zeros(1,k+1);
  122.     for i =0:k
  123.         RVector(i+1)=correlation(samp,i)
  124.     end
  125. endfunction
  126.  
  127. //получение вектора нормированных корреляционных функций от R(0) до R(k)
  128. function nRVector=normCorrelationVector(samp, k)
  129.     nRVector=zeros(1,k+1);
  130.     for i =0:k
  131.         nRVector(i+1)=normCorrelation(samp,i)
  132.     end
  133. endfunction
  134.  
  135.  
  136. //2. Построение моделей АР(M) и выбор наилучшей модели
  137. //Нахождение коэффициентов betta и alpha
  138. function betas=foundBetas(M, N)
  139.     betas = zeros(1, M);
  140.     A = zeros(M, M);
  141.     b = zeros(1, M);
  142.     for i = 1 : M
  143.         for j = 1 : M
  144.             A(i, j) = correlation(x,N + i - j);
  145.         end
  146.         b(i) = -correlation(x,N + i);
  147.     end
  148.     betas = linsolve(A, b');
  149. endfunction
  150.  
  151. //Нахождение коэффициента alpha
  152. function alpha=calcAlpha(betas)
  153.         alpha = correlation(x,0);
  154.         M=size(betas,1);
  155.         for i = 1 : M
  156.             alpha = alpha - betas(i) * correlation(x,i);
  157.         end
  158.         alpha = sqrt(alpha);
  159. endfunction
  160.  
  161. //проверка построенной модели на устойчивость
  162. function result=checkStability(betas)
  163.     sizeB = length(betas);
  164.     if (sizeB == 0) then result = %T;
  165.     elseif ((sizeB == 1) & (abs(betas(1)) < 1)) then result = %T;
  166.     elseif ((sizeB == 2) & (abs(betas(2)) < 1) & (abs(betas(1)) < 1 - betas(2))) then result = %T;
  167.     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;
  168.     else
  169.         result = %F;
  170.     end
  171. endfunction
  172.  
  173. //вычисление теоретической корреляционной функции R(k)
  174. function tR=theorCorrelation(RX, alphas, betas, k)
  175.     nm = length(alphas)+length(betas)-1;
  176.     k = abs(k);
  177.     if (k > nm) then
  178.       tR = 0;
  179.       M = length(betas);
  180.       for j = 1 : M,
  181.         tR = tR + betas(j) * theorCorrelation(RX,alphas,betas, k - j);
  182.       end;
  183.     else
  184.       tR = RX(k+1);
  185.     end
  186. endfunction;
  187.  
  188. //Вычисление нормальной теоретической корреляционной функции r(k)
  189. //(используется так же для пунктов 3,4,5,6)
  190. function tnR=theorNormCorrelation(RX, alphas, betas, k)
  191.     tnR=zeros(k+1);
  192.     for i = 0:k
  193.         tnR(i+1) = theorCorrelation(RX,alphas,betas,i);
  194.     end;
  195.     R=theorCorrelation(RX,alphas,betas,0);
  196.     tnR = tnR./R;
  197. endfunction
  198.  
  199. //вычисление среднего квадратичного отклонения
  200. //(используется так же для пунктов 3,4,5,6)
  201. function epsilon=epsilonForNCF(r, tr)
  202.     epsilon = 0;
  203.     m = min(length(r), length(tr));
  204.     for j = 1 : m,
  205.         epsilon = epsilon + (r(j) - tr(j))^2;
  206.     end;
  207. endfunction;
  208.  
  209. //Вывод на экран всех расчётов
  210. function minM=Program2(M)
  211.     printf("ЗАДАНИЕ 2------------------------------\n");
  212.     qE=zeros(M);
  213.     for i =1:M
  214.         betas = foundBetas(i, 0);
  215.         alpha=calcAlpha(betas)
  216.         printf("Коэффициенты betta для M=" + INT_FORMAT + ":\n", i);
  217.         if (checkStability(betas)) then printf(FLOAT_FORMAT + "\n", betas);
  218.             else printf("Система нестабильна");
  219.         end
  220.         printf("Коэффициент alpha для M=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, alpha);
  221.         t=normCorrelationVector(x,10);
  222.         tr=theorNormCorrelation(correlationVector(x,10),alpha, betas, 10);
  223.         qE(i)=epsilonForNCF(t,tr);
  224.         printf("Среднее квадратичное отклонение для M=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, qE(i));
  225.     end
  226.     minM=qE(1);
  227.     for i=2:M
  228.         if minM>qE(i) then
  229.             minM=qE(i);
  230.             minM=i;
  231.         end
  232.     end
  233.     printf("Лучшей моделью является модель АР(" + INT_FORMAT + ")\n", minM);
  234. endfunction
  235.  
  236.  
  237. //3. Построение моделей СС(N) и выбор наилучшей модели
  238. //создание нелинейной системы для нахождения коэффициентов alpha
  239. function y=createNoLinearSystem(alphas, N, R)
  240.   y = zeros(N + 1, 1);
  241.   for i = 1 : N + 1
  242.     for j = 1 : (N - i + 2)
  243.       y(i) = y(i) + alphas(j) * alphas(j + i - 1);
  244.     end
  245.     y(i) = y(i) - R(i);
  246.   end
  247. endfunction
  248.  
  249. //Нахождение коэффициентов alpha
  250. function alphas=calcAlphas(N)
  251.     R=correlationVector(x,10);
  252.     alphas = zeros(N + 1, 1);
  253.     alphas = fsolve(alphas, list(createNoLinearSystem, N, R));
  254. endfunction
  255.  
  256. //Вывод на экран всех расчётов
  257. function minN=Program3(N)
  258.     printf("ЗАДАНИЕ 3------------------------------\n");
  259.     qE=zeros(N+1);
  260.     for i =0:N
  261.         alphas = calcAlphas(i);
  262.         printf("Коэффициенты alpha для N=" + INT_FORMAT + ":\n", i);
  263.         printf(FLOAT_FORMAT + "\n", alphas);
  264.         t=normCorrelationVector(x,10);
  265.         tr=theorNormCorrelation(correlationVector(x,10),alphas, [], 10);
  266.         qE(i+1)=epsilonForNCF(t,tr);
  267.         printf("Среднее квадратичное отклонение для N=" + INT_FORMAT + ":"+ FLOAT_FORMAT + "\n", i, qE(i+1));
  268.     end
  269.     minN=qE(1);
  270.     for i=2:N
  271.         if minN>qE(i) then
  272.             minN=qE(i);
  273.             minN=i;
  274.         end
  275.     end
  276.     printf("Лучшей моделью является модель CC(" + INT_FORMAT + ")\n", minN);
  277. endfunction
  278.  
  279. //4 задание
  280. R = fscanfMat(FILE_CORR); //вектор корреляционной функции R(10)
  281.  
  282. function f = arma11(p)
  283.     a0 = p(1); a1 = p(2); b1 = p(3);
  284.     x0 = p(4); x1 = p(5);
  285.     f(1) = b1*R(2) + a0*x0 + a1*x1 - R(1);
  286.     f(2) = b1*R(1) + a1*x0 - R(2);
  287.     f(3) = b1*R(2) - R(3);
  288.     f(4) = a0 - x0;
  289.     f(5) = b1*x0 + a1 - x1;
  290. end;
  291.  
  292. function f = arma12(p)
  293.     a0 = p(1);  a1 = p(2);  a2 = p(3);  b1 = p(4);  x0 = p(5);  x1 = p(6);  x2 = p(7);
  294.     f(1) = b1*R(2) + a0*x0 + a1*x1 + a2*x2 - R(1);
  295.     f(2) = b1*R(1) + a1*x0 + a2*x1 - R(2);
  296.     f(3) = b1*R(2) + a2*x0 - R(3);
  297.     f(4) = b1*R(3) - R(4);
  298.     f(5) = a0 - x0;
  299.     f(6) = b1*x0 + a1 - x1;
  300.     f(7) = b1*x1 + a2 - x2;
  301. end;
  302.  
  303. function f = arma13(p)
  304.     a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4);
  305.     b1 = p(5); x0 = p(6); x1 = p(7); x2 = p(8); x3 = p(9);
  306.     f(1) = b1*R(2) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
  307.     f(2) = b1*R(1) + a1*x0 + a2*x1 + a3*x2 - R(2);
  308.     f(3) = b1*R(2) + a2*x0 + a3*x1 - R(3);
  309.     f(4) = b1*R(3) + a3*x0 - R(4);
  310.     f(5) = b1*R(4) - R(5);
  311.     f(6) = a0 - x0;
  312.     f(7) = b1*x0 + a1 - x1;
  313.     f(8) = b1*x1 + a2 - x2;
  314.     f(9) = b1*x1 + a3 - x3;
  315. end;
  316.  
  317. function f = arma21(p)
  318.     a0 = p(1); a1 = p(2); b1 = p(3); b2 = p(4);
  319.     x0 = p(5); x1 = p(6);
  320.     f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 - R(1);
  321.     f(2) = b1*R(1) + b2*R(2) + a1*x0 - R(2);
  322.     f(3) = b1*R(2) + b2*R(1) - R(3);
  323.     f(4) = b1*R(3) + b2*R(2) - R(4);
  324.     f(5) = a0 - x0;
  325.     f(6) = b1*x0 + a1 - x1;
  326. end;
  327.  
  328. function f = arma22(p)
  329.     a0 = p(1); a1 = p(2); a2 = p(3); b1 = p(4);
  330.     b2 = p(5); x0 = p(6); x1 = p(7); x2 = p(8);
  331.     f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 + a2*x2 - R(1);
  332.     f(2) = b1*R(1) + b2*R(2) + a1*x0 + a2*x1 - R(2);
  333.     f(3) = b1*R(2) + b2*R(1) + a2*x0 - R(3);
  334.     f(4) = b1*R(3) + b2*R(2) - R(4);
  335.     f(5) = b1*R(4) + b2*R(3) - R(5);
  336.     f(6) = a0 - x0;
  337.     f(7) = b1*x0 + a1 - x1;
  338.     f(8) = b1*x1 + b2*x0 + a2 - x2;
  339. end;
  340.  
  341. function f = arma23(p)
  342.     a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4); b1 = p(5);
  343.     b2 = p(6); x0 = p(7); x1 = p(8); x2 = p(9); x3 = p(10);
  344.     f(1) = b1*R(2) + b2*R(3) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
  345.     f(2) = b1*R(1) + b2*R(2) + a1*x0 + a2*x1 + a3*x2 - R(2);
  346.     f(3) = b1*R(2) + b2*R(1) + a2*x0 + a3*x1 - R(3);
  347.     f(4) = b1*R(3) + b2*R(2) + a3*x0 - R(4);
  348.     f(5) = b1*R(4) + b2*R(3) - R(5);
  349.     f(6) = b1*R(5) + b2*R(4) - R(6);
  350.     f(7) = a0 - x0;
  351.     f(8) = b1*x0 + a1 - x1;
  352.     f(9) = b1*x1 + b2*x0 + a2 - x2;
  353.     f(10) = b1*x2 + b2*x1 + a3 - x3;
  354. end;
  355.  
  356. function f = arma31(p)
  357.     a0 = p(1); a1 = p(2); b1 = p(3); b2 = p(4);
  358.     b3 = p(5); x0 = p(6); x1 = p(7);
  359.     f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 - R(1);
  360.     f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 - R(2);
  361.     f(3) = b1*R(2) + b2*R(1) + b3*R(2) - R(3);
  362.     f(4) = b1*R(3) + b2*R(2) + b3*R(1) - R(4);
  363.     f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
  364.     f(6) = a0 - x0;
  365.     f(7) = b1*x0 + a1 - x1;
  366. end;
  367.  
  368. function f = arma32(p)
  369.     a0 = p(1); a1 = p(2); a2 = p(3); b1 = p(4); b2 = p(5);
  370.     b3 = p(6); x0 = p(7); x1 = p(8); x2 = p(9);
  371.     f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 + a2*x2 - R(1);
  372.     f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 + a2*x1 - R(2);
  373.     f(3) = b1*R(2) + b2*R(1) + b3*R(2) + a2*x0 - R(3);
  374.     f(4) = b1*R(3) + b2*R(2) + b3*R(1) - R(4);
  375.     f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
  376.     f(6) = b1*R(5) + b2*R(4) + b3*R(3) - R(6);
  377.     f(7) = a0 - x0;
  378.     f(8) = b1*x0 + a1 - x1;
  379.     f(9) = b1*x1 + b2*x0 + a2 - x2;
  380. end;
  381.  
  382. function f = arma33(p)
  383.     a0 = p(1); a1 = p(2); a2 = p(3); a3 = p(4); b1 = p(5); b2 = p(6);
  384.     b3 = p(7); x0 = p(8); x1 = p(9); x2 = p(10); x3 = p(11);
  385.     f(1) = b1*R(2) + b2*R(3) + b3*R(4) + a0*x0 + a1*x1 + a2*x2 + a3*x3 - R(1);
  386.     f(2) = b1*R(1) + b2*R(2) + b3*R(3) + a1*x0 + a2*x1 + a3*x2 - R(2);
  387.     f(3) = b1*R(2) + b2*R(1) + b3*R(2) + a2*x0 + a3*x1 - R(3);
  388.     f(4) = b1*R(3) + b2*R(2) + b3*R(1) + a3*x0 - R(4);
  389.     f(5) = b1*R(4) + b2*R(3) + b3*R(2) - R(5);
  390.     f(6) = b1*R(5) + b2*R(4) + b3*R(3) - R(6);
  391.     f(7) = b1*R(6) + b2*R(5) + b3*R(4) - R(7);
  392.     f(8) = a0 - x0;
  393.     f(9) = b1*x0 + a1 - x1;
  394.     f(10) = b1*x1 + b2*x0 + a2 - x2;
  395.     f(11) = b1*x2 + b2*x1 + b3*x0 + a3 - x3;
  396. end;
  397.  
  398.  
  399.  
  400. function [koef, bestM, bestN] = arma_coefficients(R)
  401.   r = R / R(1);
  402.   qEr = zeros(1, 9);
  403.   printf('\nCreate and research ARMA model:\n');
  404.   printf('%8s  %8s  %8s  %8s  %8s  %8s  %8s  ', 'a0', 'a1', 'a2', 'a3', 'b1', 'b2', 'b3');
  405.   printf('\n');
  406.  
  407.   [answer,fun,info(1)] = fsolve([sqrt(R(1));0;0;0;0], arma11);
  408.   if (info(1) <> 4 && norm(fun) < 1e-4)
  409.     koef(1, :) = [answer(1), answer(2), 0, 0, answer(3), 0, 0];
  410.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(1, :));
  411.     alphas = koef(1, 1:2)
  412.     betas = koef(1, 5);
  413.    
  414.     tr = theorNormCorrelation(R, alphas, betas, 10);
  415.     qEr(1) = epsilonForNCF(r, tr);
  416.   else
  417.     info(1) = 0;
  418.     printf('ARMA(1,1) does not exists');
  419.     koef(1, :) = zeros(1,7);
  420.     qEr(1) = %inf;    
  421.   end
  422.   printf('\n');
  423.  
  424.   [answer,fun,info(2)] = fsolve([sqrt(R(1));0;0;0;0;0;0], arma12);
  425.   if (info(2) <> 4 && norm(fun) < 1e-4)
  426.     koef(2, :) = [answer(1), answer(2), answer(3), 0, answer(4), 0, 0];
  427.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(2, :));
  428.     alphas = koef(2, 1:3)
  429.     betas = koef(2, 5);
  430.    
  431.     tr = theorNormCorrelation(R, alphas, betas, 10);
  432.     qEr(2) = epsilonForNCF(r, tr);
  433.   else
  434.     info(2) = 0;
  435.     printf('ARMA(1,2) does not exists');
  436.     koef(2, :) = zeros(1,7);
  437.     qEr(2) = %inf;
  438.   end
  439.   printf('\n');
  440.  
  441.   [answer,fun,info(3)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0], arma13);
  442.   if (info(3) <> 4 && norm(fun) < 1e-4)
  443.     koef(3, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), 0, 0];
  444.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(3, :));
  445.     alphas = koef(3, 1:4)
  446.     betas = koef(3, 5);
  447.    
  448.     tr = theorNormCorrelation(R, alphas, betas, 10);
  449.     qEr(3) = epsilonForNCF(r, tr);
  450.   else
  451.     info(3) = 0;
  452.     printf('ARMA(1,3) does not exists');
  453.     koef(3, :) = zeros(1,7);
  454.     qEr(3) = %inf;
  455.   end
  456.   printf('\n');
  457.  
  458.   [answer,fun,info(4)] = fsolve([sqrt(R(1));0;0;0;0;0], arma21);
  459.   if (info(4) <> 4 && norm(fun) < 1e-4)
  460.     koef(4, :) = [answer(1), answer(2), 0, 0, answer(3), answer(4), 0];
  461.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(4, :));
  462.     alphas = koef(4, 1:2);
  463.     betas = koef(4, 5:6);
  464.    
  465.     tr = theorNormCorrelation(R, alphas, betas, 10);
  466.     qEr(4) = epsilonForNCF(r, tr);
  467.   else
  468.     info(4) = 0;
  469.     printf('ARMA(2,1) does not exists');
  470.     koef(4, :) = zeros(1,7);
  471.     qEr(4) = %inf;
  472.   end
  473.   printf('\n');
  474.  
  475.   [answer,fun,info(5)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0], arma22);
  476.   if (info(5) <> 4 && norm(fun) < 1e-4)
  477.     koef(5, :) = [answer(1), answer(2), answer(3), 0, answer(4), answer(5), 0];
  478.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(5, :));
  479.     alphas = koef(5, 1:3);
  480.     betas = koef(5, 5:6);
  481.    
  482.     tr = theorNormCorrelation(R, alphas, betas, 10);
  483.     qEr(5) = epsilonForNCF(r, tr);
  484.   else
  485.     info(5) = 0;
  486.     printf('ARMA(2,2) does not exists');
  487.     koef(5, :) = zeros(1,7);
  488.     qEr(5) = %inf;
  489.   end
  490.   printf('\n');
  491.  
  492.   [answer,fun,info(6)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0;0], arma23);
  493.   if (info(6) <> 4 && norm(fun) < 1e-4)
  494.     koef(6, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), answer(6), 0];
  495.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(6, :));
  496.     alphas = koef(6, 1:4);
  497.     betas = koef(6, 5:6);
  498.    
  499.     tr = theorNormCorrelation(R, alphas, betas, 10);
  500.     qEr(6) = epsilonForNCF(r, tr);
  501.   else
  502.     info(6) = 0;
  503.     printf('ARMA(2,3) does not exists');
  504.     koef(6, :) = zeros(1,7);
  505.     qEr(6) = %inf;
  506.   end
  507.   printf('\n');
  508.  
  509.   [answer,fun,info(7)] = fsolve([sqrt(R(1));0;0;0;0;0;0], arma31);
  510.   if (info(7) <> 4 && norm(fun) < 1e-4)
  511.     koef(7, :) = [answer(1), answer(2), 0, 0, answer(3), answer(4), answer(5)];
  512.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(7, :));
  513.     alphas = koef(7, 1:2);
  514.     betas = koef(7, 5:7);
  515.    
  516.     tr = theorNormCorrelation(R, alphas, betas, 10);
  517.     qEr(7) = epsilonForNCF(r, tr);
  518.   else
  519.     info(7) = 0;
  520.     printf('ARMA(3,1) does not exists');
  521.     koef(7, :) = zeros(1,7);
  522.     qEr(7) = %inf;
  523.   end
  524.   printf('\n');
  525.  
  526.   [answer,fun,info(8)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0], arma32);
  527.   if (info(8) <> 4 && norm(fun) < 1e-4)
  528.     koef(8, :) = [answer(1), answer(2), answer(3), 0, answer(4), answer(5), answer(6)];
  529.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(8, :));
  530.     alphas = koef(8, 1:3);
  531.     betas = koef(8, 5:7);
  532.    
  533.     tr = theorNormCorrelation(R, alphas, betas, 10);
  534.     qEr(8) = epsilonForNCF(r, tr);
  535.   else
  536.     info(8) = 0;
  537.     printf('ARMA(3,2) does not exists');
  538.     koef(8, :) = zeros(1,7);
  539.     qEr(8) = %inf;
  540.   end
  541.   printf('\n');
  542.  
  543.   [answer,fun,info(9)] = fsolve([sqrt(R(1));0;0;0;0;0;0;0;0;0;0], arma33);
  544.   if (info(9) <> 4 && norm(fun) < 1e-4)
  545.     koef(9, :) = [answer(1), answer(2), answer(3), answer(4), answer(5), answer(6), answer(7)];
  546.     printf('%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f ', koef(9, :));
  547.     alphas = koef(9, 1:4);
  548.     betas = koef(9, 5:7);  
  549.     tr = theorNormCorrelation(R, alphas, betas, 10);  
  550.     //disp(r);
  551.     //disp(tr);
  552.     qEr(9) = epsilonForNCF(r, tr);
  553.   else
  554.     info(9) = 0;
  555.     printf('ARMA(3,3) does not exists');
  556.     koef(9, :) = zeros(1,7);
  557.     qEr(9) = %inf;
  558.   end
  559.   printf('\n');
  560.  
  561. for i=1:9
  562.     if (qEr(i) == %inf)
  563.         printf('Inf\n');
  564.     else
  565.         printf('%.4f\n', qEr(i));
  566.     end
  567. end
  568.   printf('\n');
  569.  
  570. //выберем модель с минимальным значением qEr - это будет наилучшая модель
  571.   minqEr = min(qEr);
  572.   best = -1;
  573.   for i = 1:9
  574.       if (qEr(i) == minqEr)
  575.           best = i;
  576.       end
  577.   end  
  578.   bestM = floor((best - 1)/3) + 1; //floor - функция округления вниз
  579.   bestN = modulo(best - 1, 3) + 1;
  580.   koef = koef(best,:);
  581. endfunction
  582.  
  583. [coefARMA, armaM, armaN] = arma_coefficients(R);
  584.  
  585. function res = createProcessValues(coefProcess, armaM, armaN)
  586. r=rand(1, 11000,"normal"); //генерируем 11000 значений, распределенных по нормальному закону
  587. res = zeros(1, 10000);
  588. ksi = r(1001:11000);
  589. for i = 1:10000
  590.     j = 1;
  591.     //этта
  592.     while (j <= armaM && i > j)
  593.         res(i) = res(i) + coefARMA(4+j) * res(i-j); //coefARMA(4+j); берем 4, т.к. беты записаны 5, 6, 7 элементом
  594.         j = j + 1;
  595.     end
  596.     d = 1;
  597.     //кси
  598.     while (d <= armaN && i >= d)
  599.         res(i) = res(i) + coefARMA(d) * ksi(i-d+1); //coefARMA(d); берем d сразу, т.к альфы записаны под 1, 2, 3, 4
  600.         d = d + 1;
  601.     end
  602.     //математическое ожидание
  603.     sumTemp = 0;
  604.     for k = 1:armaM
  605.         sumTemp = sumTemp + coefARMA(4+k);
  606.     end
  607.     res(i) = res(i) + mX * (1 - sumTemp);        
  608. end
  609. endfunction
  610.  
  611. resARMA = createProcessValues(coefARMA, armaM, armaN);
  612.  
  613.  
  614. //запишем коэффициенты наилучшей модели в переменную
  615. minAR = Program2(3);
  616. minMA = Program3(3);
  617. betasAR = foundBetas(minAR, 0);
  618. alphaAR=calcAlpha(betasAR);
  619. coefAR = [alphaAR(1), 0, 0, 0, betasAR(1), betasAR(2), betasAR(3)];
  620. alphaMA = calcAlphas(minMA);//для СС
  621. coefMA = [alphaMA(1), alphaMA(2), alphaMA(3), alphaMA(4), 0, 0, 0];
  622. resAR = createProcessValues(coefAR, minAR, 0);
  623. resMA = createProcessValues(coefMA, 0, minMA);
  624.  
  625. //построение графиков для 5 задания
  626. //график для АРСС
  627. function processPlotForARMA()
  628.     t = [0:10];
  629.     alpha = coefARMA(1:4);
  630.     betas = coefARMA(5:7);  
  631.     normCorr_x=normCorrelationVector(x,10);
  632.     theorNormCorr=theorNormCorrelation(correlationVector(x,10),alpha, betas, 10);
  633.     normCorr_res=normCorrelationVector(resARMA',10);
  634.     F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
  635.     plot2d(t,normCorr_x, axesflag = 1,style = 5);  
  636.     plot2d(t,theorNormCorr, axesflag = 1,style = 3);  
  637.     plot2d(t,normCorr_res, axesflag = 1,style = 2);  
  638.     xtitle("","Индекс","Нормированная корреляционная функция");
  639.     xgrid();
  640.     legend('Процесс', 'Теоретическая АРСС(3,3)', 'Смоделированный процесс');
  641. endfunction;
  642. processPlotForARMA();
  643.  
  644. //график для АР
  645. function processPlotForAR()
  646.     t = [0:10];
  647.     normCorr_x=normCorrelationVector(x,10);
  648.     theorNormCorr=theorNormCorrelation(correlationVector(x,10),alphaAR, betasAR, 10);
  649.     normCorr_res=normCorrelationVector(resAR',10);
  650.     F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
  651.     plot2d(t,normCorr_x, axesflag = 1,style = 5);  
  652.     plot2d(t,theorNormCorr, axesflag = 1,style = 3);  
  653.     plot2d(t,normCorr_res, axesflag = 1,style = 2);  
  654.     xtitle("","Индекс","Нормированная корреляционная функция");
  655.     xgrid();
  656.     legend('Процесс', 'Теоретическая АР(3)', 'Смоделированный процесс');
  657. endfunction;
  658. processPlotForAR();
  659.  
  660. //график для СС
  661. function processPlotForMA()
  662.     t = [0:10];  
  663.     normCorr_x=normCorrelationVector(x,10);
  664.     theorNormCorr=theorNormCorrelation(correlationVector(x,10),alphaMA, [], 10);
  665.     normCorr_res=normCorrelationVector(resMA',10);
  666.     F=figure('figure_name','Оценка нормированной корреляционной функции', 'background',8);
  667.     plot2d(t,normCorr_x, axesflag = 1,style = 5);  
  668.     plot2d(t,theorNormCorr, axesflag = 1,style = 3);  
  669.     plot2d(t,normCorr_res, axesflag = 1,style = 2);  
  670.     xtitle("","Индекс","Нормированная корреляционная функция");
  671.     xgrid();
  672.     legend('Процесс', 'Теоретическая СС(3)', 'Смоделированный процесс');
  673. endfunction;
  674. processPlotForMA();
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement