Generate normal random sequence that is based on the given conditions. Also calculate different stats and show correlation graphs.

normalGenerator.m

function normalGenerator()
  count = 1000;
  actualCount = count*2;
  adCount = count * 2 + 100;
  mu = 2;
  sigma = 0.5;
  
  kSize = 20;
  t = 1.964;
  
  a = 0.1;
  k = 1.73 * (count ^ (1/3));
  
  floatValues = generateFloatValues(adCount);
  [float1, float2] = divideFloatValues(floatValues);
  
  X1 = zeros(1, count);
  X2 = zeros(1, count);
  for i = 1:count
    X1(i) = withOffset(getX1(float1, float2, i), mu, sigma);
    X2(i) = withOffset(getX2(float1, float2, i), mu, sigma);
  end
  
  X = [X1, X2];
  
  mat = getMat(X)
  dis = getDis(X, mat)
  confidenceInterval = getConfidenceInterval(X, dis, mat, t)
  showCorrelation(X, mat, kSize);
  showNextPreviousCorrelation(X);
  showFrequency(X, -3, 7, 60);
  
  [Z, zTable] = getZ(mu, sigma, X, k, actualCount);
  zTable
  Z
  X2 = calcX2(a, k)
end

function [Z, zTable] = getZ(mu, sigma, X, k, actualCount)
  p = zeros(k, 1);
  h = zeros(k, 1);
  normal = @(x) exp(-((x-mu).^2)/(2*sigma.^2))/sqrt(2*pi*sigma.^2);
  delta = (max(X)-min(X))/k;
  from = min(X);
  to = from + delta;
  for i = 1:k
    p(i) = quad(normal, from, to);
    h(i) = countBetween(from, to, X);
    from = from + delta;
    to = to + delta;
  end
  
  Z = 0;
  zTable = zeros(k, 3);
  for i = 1:k
    localZ = ((h(i) - actualCount*p(i))^2)/(actualCount*p(i));
    Z = Z + localZ;
    
    zTable(i, 1) = localZ;
    zTable(i, 2) = h(i);
    zTable(i, 3) = p(i);
  end
end

function X2 = calcX2(a, k)
  X2 = chi2inv(1-a, k-1);
end

function res = countBetween(from, to, X)
  [xSize, ySize] = size(X);
  res = 0;
  for j = 1:ySize
    if (X(j) > from && X(j) < to)
      res = res+1;
    end
  end
end
    
function discrete = getDiscreteFromFloat(value, inputData, intervals)
  for i = 1:size(intervals)
    if value < intervals(i)
      discrete = inputData(1, i);
      break;
    end
  end
end

function xWithOffset = withOffset(x, mu, sigma)
  xWithOffset = mu + sigma*x;
end

function x = getX1(y1, y2, i)
  x = sqrt(-2*log(y1(i)))*sin(2*pi*y2(i));
end

function x = getX2(y1, y2, i)
  x = sqrt(-2*log(y1(i)))*cos(2*pi*y2(i));
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
  y = y(2:count);
end