lab3

function MC_LR_3_KHMELEV
  V = 17
  a_ = -0.7; b_ = 7.3; % по условию


 %_______________________________ЗАДАНИЕ I______________________________________

  xx = [];
  xx = sort(oneOfArrays('nd')); % функция, выдающая выборку в соответствии с вариантом, в самом низу
  [N, m] = SERIES(1, xx);
  printf('\nThe number of items in our array: N = %i\nThe number of intervals: m = %i\n\n', N, m);

  Table(1); % первая таблица


  % график плотности нормального распределения, наложенный на гистограмму относительных частот:
  for j = 1:100
    b(j) = a(1) + j*d/100;
    c(j) = normpdf((b(j)-A)/sigma)/sigma;
  endfor
  figure('Name','Normal probability density function & histogram');
  bar(x, w/h, 'hist'); % график гистограммы относительных частот
   hold on % задерживаем предыдущий график, чтобы наложить на него график плотности:
  plot(b, c, "r");
   grid on % включаем сеточку
   hold off


  Table(2); % вторая таблица


  % проверка гипотезы о соответствии выборки нормальному распределению:
  X = myround(sum(NWP), 5); % наше значение критерия x
  Hypothesis(1, X);





 %_______________________________ЗАДАНИЕ II_____________________________________

  xx = [];
  xx = sort(oneOfArrays('ud'));
  [N, m] = SERIES(2, xx);
  printf('\nThe number of items in our array: N = %i\nThe number of intervals: m = %i\n\n', N, m);


  % график плотности равномерного распределения на отрезке [a, b], наложенный на гистограмму относительных частот:
  figure('Name','Uniform probability density function & histogram');
  bar(x, w/h, 'hist');
   hold on % задерживаем предыдущий график, чтобы наложить на него график плотности:
  plot([x(1)-h/2, x(end)+h/2], [1/(b_-a_), 1/(b_-a_)], 'r');
   grid on, hold off


  Table(3); % третья таблица


  % проверка гипотезы о соответствии выборки равномерному распределению на отрезке [a, b]:
  X = myround(sum(NWP), 5);
  Hypothesis(2, X);





 %______________________________________________________________________________

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

  function [N, m] = SERIES(task, xx) % ряд и гистограмма
    m = N = d = h = n = k = 0;
    a = w = x = [];
    printf('\n\n-------------\n  TASK NUMBER %i\n\n\n', task);

    N = length(xx);
    m = 1 + fix(log2(N)); % число интервалов по формуле Стерджеса (fix отбрасывает дробную часть)

    if task == 1
      a(1) = xx(1); %а0
      a(m+1) = xx(end); %аm
    else
      a(1) = a_; %а0
      a(m+1) = b_; %аm
    endif

    d = a(m+1) - a(1);
    h = d/m; % шаг разбиения
    for k = 1:m-1
      a(k+1) = a(k) + h; %от а0 до аm (но у нас от а1 до а(m+1))
    endfor

    [n, x] = hist(xx, m); % находим число значений xx, попавших в интервал [аk, аk+1] - это n[]; и середины таких интервалов - это x[]
    w = n/N; % значения относительных частот

    printf('Interval series:\n\n['); % группированная выборка (интервальный вариационный ряд);
    for i = 1:m
      printf('%i, %i]\n n = %i, w = %i\n\n', a(i), a(i+1), n(i), myround(w(i), 5));
      if i == m
        printf('\n');
        break
      endif
      printf('(');
    endfor


    A = 0; % мат. ожидание
    A = mathExpect(w, x, a, 1);

    SIGMA2 = 0; % дисперсия
    SIGMA2 = mathExpect(w, x, a, 2) - A^2 - h^2/12;

    sigma = 0;
    sigma = sqrt(SIGMA2);

    if task == 1 % т.к. вывести эти значения требуется только в задании I
      printf('Approx. Math. expectation: a = %i\nApprox. Variance: sigma^2 = %i\n', A, SIGMA2);
    endif

  endfunction


  function Table(number) % построение таблиц 1, 2 (для задания I) и 3 (для задания II)
    p = [];
    printf('\n------\n\n Table %i:\n\n', number);

    if number == 3 % для таблицы 3:
      p = ones(1, m+1); % p(1) (или p(0), как указывается в задании) мы все равно потом выкинем (т.к. в задании счет начинается с 1, а не с 0), и их тогда будет m, а не m+1
      p *= 1/m;
    else % для таблиц 1 и 2:
      t = (a - A)/sigma;
      p(2) = Fi(t(2));
      p(3:m) = Fi(t(3:m)) - Fi(t(2:m-1));
      p(m+1) = 1 - Fi(t(m)); % та же ситуация, что и описанная чуть выше
    endif


    if number == 1 % таблица 1:
      TABLE = [];

      f = fi(t)/sigma;
      F = Fi(t);

      for k = 1:m+1 % индексы выводятся в соответствии с примером таблицы в задании
        TABLE(k, :) = [k-1, a(k), t(k), f(k), F(k), p(k)];
      endfor

      printf('      k        ak    (ak-A)/sgm  fi/sgm      Fi        pk\n');
      disp(TABLE);

    else % таблицы 2 и 3:
      WP = NWP = p(1) = []; % выкидываем p(1) (который по заданию является p(0))

      WP = abs(w - p);
      NWP = N.*WP.^2 ./ p; % поэлементно возводим в квадрат и так же поэлементно делим на элементы массива p

      printf('   k   [a(k-1), a(k)]        wk       pk    |wk-pk|   N(wk-pk)^2/p\n\n');
      for k = 1:m % индексы будут выводиться такими, какие требуется в задании (то есть без учета того, что в программе мы везде начинаем с единицы, а не с нуля)
        printf('   %i [%i, %i]   %i   %i   %i     %i\n', k, a(k), myround(a(k+1), 5), myround(w(k), 5), myround(p(k), 5), myround(WP(k), 5), myround(NWP(k), 5));
      endfor
      printf('\n max|wk-pk| = %i\n sum(N(wk-pk)^2/p) = %i\n', myround(max(WP), 5), myround(sum(NWP), 5));
      printf('\n-------\n\n');
    endif

  endfunction


  function Hypothesis(task, X) % проверка гипотезы о соответствии выборки
    aa = 0.05; % aa - уровень значимости
    Xkp(4:8) = [9.5, 11.1, 12.6, 14.1, 15.5] % таблица критических значений Xкр

    if task == 1
      l = m - 3; % l - число степеней свободы для задания I
    else
      l = m - 1; % для задания II
    endif

    if X <= Xkp(l)
      printf('X = %i; Xkp(m-3) = %i, %i<=%i  =>  Gipoteza mozhet byt prinyata\n', X, Xkp(l), X, Xkp(l));
    else
      printf('X = %i; Xkp(m-3) = %i, %i>%i  =>  Gipoteza nie mozhet byt prinyata\n', X, Xkp(l), X, Xkp(l));
    endif

  endfunction

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

end




%_______________________________________________________________________________

%%%%%%%%%%%%%%%%%%% ОСТАЛЬНЫЕ ФУНКЦИИ %%%%%%%%%%%%%%%%%%%


function A = mathExpect(w, x, a, t); %математическое ожидание
    A = 0;

    A = sum(w.*x.^t); % x - это x* из условия МС_ЛР_3, уже посчитанный функцией hist()

endfunction


function fi = fi(t) % плотность стандартного нормального распределения N(0,1)
  fi = 0;  %f = 1/sqrt(pi = 2*3.1415926535)*exp(-t^2/2)

  fi = normpdf(t, 0, 1);

endfunction

function Fi = Fi(t) % функция распределения стандартного нормального закона N(0,1)
  Fi = 0;

  Fi = normcdf(t, 0, 1);

endfunction


%function A = mathExpect(w, m, a, t); %математическое ожидание, посчитанное без встроенных функций)
%   A = 0;
%
%   for k = 1:m
%       xk = ((a(k) + a(k+1))/2)^t;
%       A += xk*w(k);
%  endfor
%
%endfunction


function R = myround(x, k)
  R = 0;

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

endfunction


function arr = oneOfArrays(a)
  arr = []; %вариант 17

    if a == 'nd'
        arr = [ 2.34077, -0.47612, -0.59152, -1.33363, -0.63093, 0.48788, -4.95074, 0.56767, 1.19521, -2.79534, -1.20603, 1.30688, -2.85712, -0.55494, -1.90109, 3.60836, 0.21293, -1.17190, -1.21709, -0.86380, 2.78199, -1.12639, 0.16847, -0.15564, -0.04350, -3.22104, -1.94865, -1.30119, 0.14260, -0.97421, -1.37667, -7.28036, 4.31504, -2.30376, -4.21841, -4.12702, -2.45313, -0.68081, 0.65864, -0.38498, -2.20725, 1.80698, -1.46469, -2.83590, -1.91566, 2.85141, -3.82743, 3.16380, -1.75219, 1.81357, -2.26096, -3.19839, -0.37549, -1.86859, -0.31260, -1.28324, 2.38648, 1.75209, -1.12582, 3.20351, -0.85097, 0.61758, 2.75101, -1.15348, -3.24758, -0.62437, -0.31909, 0.95708, -1.08075, -1.05882, 0.32362, -5.49754, -3.43182, 0.95152, 2.93082, -0.44737, -0.60788, 0.25537, -2.80409, -2.22244, -6.27760, 0.76520, 3.15344, 1.98039, -4.35556, -2.92046, 3.17636, -1.54461, 3.19059, -0.62366, -2.41449, 2.34206, -0.56942, -1.25159, -2.00228, -0.72641, -3.44084, -2.72999, 1.97637, -5.51718, -0.04327, -1.63052, -1.65610, -1.57951, -4.70447, 2.95226, 1.90369, -2.50510, 2.34484, -1.85043, 1.79099, -1.45891, 3.46837, 1.89493, -2.40129, -3.11331, -0.95251, 0.99113, 1.88399, -0.42847, -2.58520, 1.38002, 1.22401, -2.92613, -0.07367, 1.31508, 2.17188, 0.24897, -0.70232, -0.16624, -1.26202, -3.24450, -1.17668, 0.84247, -0.59938, -2.10478, 1.74759, -0.00428, -1.54944, 0.16808, -0.12737, -3.62850, 1.72206, 0.69783, 1.25305, -2.32759, -1.58966, -2.53622, 4.07967, -1.01046, -0.29265, 0.08308, -0.91494, 4.14844, -2.13790, 4.27977, 0.97369, -3.21341, -1.90497, 1.57068, 0.90532, 1.32906, 0.00284, 0.30189, 1.55175, 0.50187, 2.82884, -6.23639, 0.48637, 2.99302, 2.01778, -4.05521, -3.84553, -1.39976, -2.65657, -0.57231, -5.44073, -0.01255, 2.45312, 0.38175, 0.33639, 0.10010, 0.56131, -0.94303, -5.03709, 0.31000, -2.87028, 0.51529, -1.88752, -0.48789, -5.28494, -0.91663, -0.05016, -0.60857, -1.48554, -0.71420, -1.29857, 1.67761, -2.78399, 0.55937, -0.48658, -0.16625, -1.35844, -2.03059, 2.35581, -0.40343, 1.90275, -1.30859, -0.74700, -3.78552, -4.86740, -0.98373, 0.21135, -0.05607, -4.51215, -0.37222, 2.70644, -1.46749, -0.36272, -0.28096, 2.05548, 1.03205, -3.97384, -3.26241, -1.67013, -4.65148, 1.31498, -2.49333, -1.70760, -1.31645, 2.15126, 1.01024, 1.36154, 0.70141, -3.75326, -4.25257, 0.61316, -4.97811, -0.64413, 0.54008, -0.15627, 0.95668, 2.33636, -0.10633, 0.18355, -0.66052, -2.62649, -0.33295, 1.16275, -4.06517, -3.07452, 0.03409, -0.53623, -2.65940, 3.00544, -0.52499, 0.59608, 3.77866, -2.25218, 1.91575, 1.46961, -1.45929, -0.89563, -0.85259, -3.96968, -0.97911, 2.14339, -5.94468, 2.08057, 0.47641 ];
    elseif a == 'ud'
        arr = [ 5.16392, 6.12984, 2.13544, 1.32408, 1.96856, 4.02600, 5.75448, 6.16088, 3.84952, 1.29904, 5.16816, 5.95720, 2.96240, 4.46312, 0.09080, 3.13536, 2.61328, 3.72544, 3.12392, 5.96696, 1.70104, 3.23016, 1.38256, 2.54120, 4.25520, 6.30704, 1.66768, 2.31632, 0.20840, 0.20296, 2.58312, 1.08168, 2.45432, 3.51280, 0.89648, 3.02304, 6.37296, 5.66424, 1.39224, 4.70040, 5.09552, 6.04120, 5.85832, 5.41784, 3.79408, 1.73416, 6.63960, 3.81488, 0.62632, 0.83240, 6.22136, 4.78584, 1.07128, 5.59656, 2.37768, 2.48632, 3.46456, 3.48200, 1.43648, 2.90848, 2.34792, -0.52328, 3.49880, -0.34240, 5.94280, 3.56072, -0.05208, 1.12792, 3.05880, 6.83288, 4.14928, 1.81496, 1.96480, 5.04352, -0.33960, 7.05272, 3.49360, 0.47408, 6.85448, 7.12656, 5.48552, 3.92768, 3.10816, 0.91832, 3.38720, 6.47200, 5.59824, -0.52552, 2.56408, 2.25192, 0.08120, 4.83488, 2.00104, 3.81152, 0.66064, 4.36624, 2.04064, 2.63544, 0.56288, 0.83376, 0.78960, 5.15960, 0.02120, 6.17336, 4.39728, 1.18256, 6.66152, 5.04624, 7.19472, 0.44624, 2.24224, 2.19896, -0.55984, 1.10120, -0.14064, 4.36096, 6.74456, 2.94256, 5.45328, 3.12000, 4.33264, -0.69456, 4.62424, 0.68104, 0.82648, 7.10696, 6.04744, -0.31392, 3.30640, 6.71200, 1.74944, 2.65432, 6.57728, 3.64768, 6.21944, 5.12928, -0.51664, 1.76864, 3.10448, -0.12808, 5.13712, 1.92320, 4.16928, 5.75504, 0.49512, 3.92760, 1.13736, 1.55192, 3.40408, 2.81448, 5.87832, 3.92568, -0.69184, 2.42792, 4.41704, 2.25088, 1.34296, 5.88008, 4.49728, 3.83904, 1.94536, 4.41712, 5.46872, 6.07192, 3.49144, 1.26456, 4.93832, 5.22912, 5.94144, 3.07192, 0.93896, -0.35784, 5.59520, 0.85528, 3.78480, 2.85904, 4.76072, 5.17104, 0.08912, 1.15784, 6.04888, 1.98080, 6.12296, 1.32280, 6.75840, 5.90520, 3.05696, 6.50248, 4.91984, 4.93072, -0.21440, 6.97032, 7.13704, 1.19496, 3.58536, 6.22376, -0.13648, 1.02504, 0.85336, 2.32040, 6.71400, 0.99768, 2.72640, 7.20304, -0.39520, 3.55688, 0.38304, 5.81256, 1.94352, 6.03600, 5.33640, 1.34592, 1.98472, 4.61624, 3.35552, 0.22360, 6.59256, 5.38112, -0.17688, 4.79592, 4.31848, 4.55432, 3.88256, -0.25816, 5.66152, 4.66528, 6.25184, 0.69288, 6.09736, -0.64000, 1.54936, 0.50144, 4.47960, 2.77216, 6.71168, 6.76128, 6.20568, 0.74952, 4.57984, 4.05544, 0.47800, 5.18040, 4.18304, 5.86688, 6.99248, 2.74312, 2.00624, 2.58472, 0.70040, 6.15032, 2.28496, 2.77224, 0.86704, -0.46992, 3.38904, 3.39376, 2.38592, 0.34264, 5.96712, 1.13424, 1.45568, 0.87744, -0.18328, 4.22968, 3.52856, -0.05936, 3.99336, 2.05608, 2.60672, 6.97056 ];
    endif

endfunction