lab5

function MC_LR_5_KHMELEV

  V = 17
  U = oneOfArrays;
  [N, m] = size(U);
  printArr5(N, m);

 %_______________________________ЗАДАНИЕ I______________________________________

  printf('\n\n-------------\n  TASK NUMBER 1\n\n');


  % пункт 1
  printf('\n  Paragraph 1:\n\n');
  TABLE_1 = TABLE_4 = TABLE_5 = TABLE2_1 = TABLE2_3 = [];
  hyp1_1 = [];
  for i = 1:3
    % для первого пункта задания I:
    Calculations( U(:, i-mod(1,i)), U(:, 0.5*(i-mod(i,2))+2) ) %столбцы (1,2), (1,3), (2,3)
    TABLE_1(i, :) = [x, y, x2, y2, Sx, Sy, Tnn];
    % для первого пункта задания II:
    TABLE2_1(i, :) = [S1, S2, N-1, N-1, Fnn];
    % для четвертого пункта задания I:
    TABLE_4 = [TABLE_4, t_test_2(U(:, i-mod(1, i)), U(:, 0.5*(i-mod(i,2))+2))]; % добавляем к предыдущим элементам массива наовые: результаты теста для U(:, 1) и (U:, 2); для U(:, 1) и (U:,3); и затем для U(:, 2) и U(:, 3)
    % для пятого пункта задания I:
    TABLE_5 = [TABLE_5, welch_test(U(:, i-mod(1, i)), U(:, 0.5*(i-mod(i,2))+2))];
    % для третьего пункта задания II:
    TABLE2_3 = [TABLE2_3, var_test(U(:, i-mod(1, i)), U(:, 0.5*(i-mod(i,2))+2))];
  endfor
  printf('       x          y         x2         y2         Sx         Sy        Tnn\n');
  disp(myround(TABLE_1, 5));
  hypothesis(1.1, TABLE_1(:, end));

  % пункт 2
  printf('\n\n  Paragraph 2:\n\n');
  S_calculation(U); % гипотеза проверяется в этой же функции

  % пункт 3
  printf('\n\n  Paragraph 3:\n\n');
  pval = anova(U)
  hypothesis(1.3, pval);

  % пункт 4
  printf('\n\n  Paragraph 4:\n\n pval[t_test_2] \n');
  disp(myround(TABLE_4', 5));
  hypothesis(1.4, TABLE_4');

  % пункт 5
  printf('\n\n  Paragraph 5:\n\n pval[welch_test] \n');
  disp(myround(TABLE_5', 5));
  hypothesis(1.5, TABLE_5');



 %_______________________________ЗАДАНИЕ II_____________________________________


  printf('\n\n-------------\n  TASK NUMBER 2\n\n');
  % все таблицы подсчитываются в цикле в задании I

  % пункт 1
  printf('\n\n  Paragraph 1:\n\n');
  printf('      S1        S2       k1        k2        Fnn\n');
  disp(myround(TABLE2_1, 5));
  hypothesis(2.1, TABLE2_1(:, 3:5));

  % пункт 2
  printf('\n\n  Paragraph 2:\n\n');
  pval = bartlett_test(U(:, 1), U(:, 2), U(:, 3))
  hypothesis(2.2, pval);

  % пункт 3
  printf('\n\n  Paragraph 3:\n\n pval[var_test] \n');
  disp(myround(TABLE2_3', 5));
  hypothesis(2.3, TABLE2_3');






 %______________________________________________________________________________

 %%%%%%% ФУНКЦИИ, ВКЛЮЧЕННЫЕ В ОСНОВНУЮ ФУНКЦИЮ %%%%%%%

  function Calculations(X, Y)

    % для задания I:
    x = summation(X, 1)/N;
    x2 = summation(X, 2)/N;
    y = summation(Y, 1)/N;
    y2 = summation(Y, 2)/N;

    Sx = N/(N - 1)*(x2 - x^2);
    Sy = N/(N - 1)*(y2 - y^2);

    Tnn = (x - y)/sqrt((Sx + Sy)*(N - 1))*sqrt(N*(N - 1)); % у нас M = N, поэтому функция выглядит проще

    % для задания II:
    S1 = max(Sx, Sy);
    S2 = min(Sx, Sy);

    Fnn = S1/S2;

  endfunction


  function S_calculation(U)

    sum_U = sum(U); % три элемента, каждый из которых это сумма столбца от i = 1 до N
    u = sum(sum_U)/(N*m); % u с чертой
    uj = sum_U/N; % массив uj с чертой

    S_total = sum(summation(U - u, 2))
    S_fact = N*summation(uj - u, 2)
    S_ost = S_total - S_fact

    k1 = m-1;
    k2 = m*(N-1);

    S2_fact = S_fact/k1
    S2_ost = S_ost/k2

    k1
    k2

    Fnm = S2_fact/S2_ost

    hypothesis(1.2, [k1, k2, Fnm]);

  endfunction


  function hypothesis(number, param)
    tf = 0;
    aa = 0.05;
    printf('\n  Proverka gipotezy (0 - NET, 1 - DA)\n\n');

    switch number

      case 1.1
        TNN = abs(param);
        tkp = tinv(0.975, 2*N-2)*ones(3, 1);

        tf = TNN <= tkp;

        printf('    |Tnn|     tkp(2N-2)    vyvod\n');
        disp(myround([TNN, tkp, tf], 5));

      case 1.2
        k1 = param(1);
        k2 = param(2);
        Fnm = param(3);
        z = finv(0.95, k1, k2)

        tf = Fnm <= z;

        printf('vyvod: %i\n', tf);

      case 1.3
        pval = param;

        tf = pval >= aa;

        printf('vyvod: %i\n', tf);

      case 1.4
        PVAL = param;
        AA = aa*ones(3, 1);

        tf = PVAL >= AA;
        printf(' pval[t_test_2]   a        vyvod\n');
        disp(myround([PVAL, AA, tf], 5));

      case 1.5
        PVAL = param;
        AA = aa*ones(3, 1);

        tf = PVAL >= AA;
        printf(' pval[welch_test]  a        vyvod\n');
        disp(myround([PVAL, AA, tf], 5));

      case 2.1
        k1 = param(:, 1);
        k2 = param(:, 2);
        FNN = param(:, 3);
        Fa = finv(0.975, k1, k2);

        tf = FNN <= Fa;
        printf('    |Fnn|     F      vyvod\n');
        disp(myround([FNN, Fa, tf], 5));

      case 2.2
        pval = param;

        tf = pval >= aa;

        printf('vyvod: %i\n', tf);

      case 2.3
        PVAL = param;
        AA = aa*ones(3, 1);

        tf = PVAL >= AA;
        printf(' pval[var_test]   a        vyvod\n');
        disp(myround([PVAL, AA, tf], 5));

    endswitch

  endfunction


  function S = summation(t, k)
    S = 0;% Суммирование элементов в степени

    S = sum(t.^k);

  endfunction

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end


function R = myround(x, k) % округление до 5 знаков после запятой
  R = 0;

  R = round(x*10^k)/10^k;

endfunction


function arr = oneOfArrays % выдача массива
  arr = []; %вариант 17

    arr = [ -2.42349, -2.03716, -1.93880;
          -3.72517, 2.81466, -1.57690;
          -3.74582, 0.74513, 3.97149;
          1.40175, -0.80514, 0.57361;
          -1.17462, -2.31128, 2.96124;
          0.32442, -1.77895, 1.09566;
          0.11940, 3.40010, -3.26145;
          0.62900, -0.39564, -2.56283;
          0.34779, -0.56814, 2.64745;
          0.14588, 3.11038, -0.80312;
          -0.41323, 1.45303, 0.38234;
          -1.34367, 0.24958, -0.51843;
          -2.80189, 2.21146, -0.92177;
          0.18032, -0.60075, 1.43413;
          1.73334, -1.12983, 2.33078;
          -0.47590, -0.61128, -1.09554;
          0.00733, 2.11087, 3.05404;
          -2.26858, 0.24760, -1.85628;
          2.08119, 1.19782, -0.50469;
          0.07452, -2.51749, 0.87703 ];

endfunction


function printArr5(N, m)
  arr = oneOfArrays;
  len = N*m;

  printf('U = \n\n');
  for i = 1:len
    printf('%i', arr(i));
    if mod(i, 3) == 0 && i ~= len
      printf(';\n');
    elseif i ~= len
      printf(', ');
    endif
  endfor

%  for j = 1:m
%    printf('\n');
%    for i = 1:N
%      printf('%i\n', arr(i, j));
%    endfor
%  endfor

endfunction