Generate Monte Karlo random sequence that is based on the given function. Also calculate different stats and show correlation graphs.

monteKarloGenerator.m

function monteKarloGenerator()
  count = 5000;
  a = 0;
  b = 4;
  M = 0.5;
  p = @func;
  
  kSize = 20;
  t = 1.964;
  %part1 = @(x) (x/4), 0, 2;
  %part2 = @(x) 0.25, 2, 4;
  %p = getP(part1, part2);
  %fplot(@(x) p(x));
  
  floatValues = generateFloatValues(count);
  [float1, float2] = divideFloatValues(floatValues);
  vals = getFloatValuesInInterval(a, b, M, float1, float2, p);
  
  mat = getMat(vals)
  dis = getDis(vals, mat)
  confidenceInterval = getConfidenceInterval(vals, dis, mat, t)
  showCorrelation(vals, mat, kSize);
  showNextPreviousCorrelation(vals);
  showFrequency(vals, a, b, 30);
end

%function p = getP(part1, part2)
%  p = @(x) part1(x)+part2(x);
%end

function y = func(x)
  if x >= 0 && x < 2
    y = x/4;
  elseif x >= 2 && x <=4
    y = 0.25;
  else
    y = 0;
  endif
end

function mat = getMat(y)
  [sizeX, sizeY] = size(y);
  mat = 0;
  for i = 1:sizeY
    mat = mat + y(i);
  end
  mat = mat/sizeY;
end

function dis = getDis(y, mat)
  [sizeX, sizeY] = size(y);
  dis = 0; 
  for i = 1:sizeY
    dis = dis + (y(i)-mat)^2; 
  end 
  dis = dis/sizeY;
end

function interval = getConfidenceInterval(y, dis, mat, t)
  [sizeX, sizeY] = size(y);
  a = mat - t*sqrt(dis/sizeY);
  b = mat + t*sqrt(dis/sizeY);
  interval = [a, b];
end

function showCorrelation(y, mat, kSize)
  [sizeX, sizeY] = size(y);
  k = zeros(kSize, 1);
  for i = 1:kSize
    for j = 1:(sizeY-i-1)
      k(i) = k(i) + (y(j)-mat)*(y(j+i-1)-mat);
    end
    k(i) = k(i)/(sizeY-i-1);
  end 
  subplot(2,2,3);
  plot(0:kSize-1,k);
  title('Correlation');
end

function showNextPreviousCorrelation(y)
  [sizeX, sizeY] = size(y);
  X = zeros(sizeY, 1);
  Y = zeros(sizeY, 1);
  for n = 1:sizeY-1
    X(n) = y(n);
    Y(n) = y(n+1);
  end
  subplot(2,2,2);
  scatter(X, Y);
  title('(y(i+1), y(i))');
end

function showFrequency(y, a, b, ints)
  [sizeX, sizeY] = size(y);
  step = (b-a)/ints;
  borders = a:step:b;
  hist_arr = zeros(ints+1, 1);
  i = 0;
  for b = borders
    q = b + step;
    i = i + 1;
    for j = 1:sizeY
      if b<y(j) && y(j)<q 
        hist_arr(i) = hist_arr(i) + 1; 
      end
    end
    hist_arr(i) = hist_arr(i)/sizeY/step;
  end   
  subplot(2,2,1);
  bar(borders, hist_arr); 
  title('histogram');
end

function vals = getFloatValuesInInterval(a, b, M, float1, float2, p)
  index = 1;
  vals = [];
  for i = 1:min(size(float1), size(float2))
    n1 = a + float1(i)*(b-a);
    n2 = M * float2(i);
    if (n2 < p(n1))
      vals(index) = n1;
      index++;
    endif
  end
end

function [float1, float2] = divideFloatValues(floatValues)
  [x, y] = size(floatValues);
  mid = round(x/2);
  float1 = floatValues(1:mid);
  float2 = floatValues(mid:x);
end

function y = generateFloatValues(count) 
  m = 2^31 - 1;
  a = 630360016;
  n = count+1;
  y = zeros(count, 1); 
  ec = zeros(count, 1);
  ec(1) = 34238443;
  
  for i = 2:1:n
      ec(i) = mod(a*ec(i-1),m);
      y(i) = ec(i)/m;
  end
end