clear;
close all;
clc;

% the clt tells you what is the distribution of the average of n random
% variables
% we will vary n from 5 to some large value to demonstrate this
% Note: make sure that you have enough ram on your machine to store a 6000
% x 10000 array!
for n=[5 10 15 50 100 200 500 1000 10000]
    nsamp = 6000; % to plot a histogram we need several samples

    
    % our sampling distribution is P(X=0)=0.35, P(X=1)=0.35, P(X=2)=0.3
    X = rand(nsamp,n); X(X < 0.35) = 0; X(X >= 0.7) = 2; X(X >= 0.35 & X < 0.7) = 1; 
    mu = 0.95; % is the true mean of this distribution (you can check this)
    sig2 = 0.6475; % the true variance of this distribution - check this yourself
    
    % what is the meaning of this nsamp?
    % suppose we generate 'n' samples from the sampling distribution - we want
    % to examine the distribution of the average value of these samples
    % for that we need to build a histogram of samples of the average value
    % itself! For this we must repeat the process of generating 'n' samples
    % followed by computation of the average value, some 'nsamp' number of
    % times!
    
    % this is what is denoted as Y_j in the slides - except that we aren't
    % dividing by sqrt(sig2) there. the distribution of the quantity below
    % should converge to that of the standard normal.
    % in the lecture slides, it should converge to N(0,sig2)
    sumX = (sum(X,2)-n*mu)/(sqrt(sig2*n));
    
    % we use 50 bins for our histogram
    numbins = 50;
   
    hist(sumX,numbins); %minval+stepsize/2:stepsize:maxval-stepsize/2);
    title (sprintf('PDF of the average of %d iid random variables',n));
    pause;
end

% we will repeat the earlier exercise with a [0,1] uniform random
% distribution - you will notice the resemblance to the standard normal
% distribution kicks in much earlier! this is generally true for 
% distributions that are symmetric about their mean.
for n=[5 10 15 50 100 200 500 1000 10000]
    nsamp = 6000; % to plot a histogram we need several samples
        
    % our sampling distribution is [0,1] uniform random
    X = rand(nsamp,n); 
    mu = 0.5; % is the true mean of this distribution (you can check this)
    sig2 = 0.0833; % the true variance of this distribution - check this yourself
    
    % this is what is denoted as Y_j in the slides - except that we aren't
    % dividing by sqrt(sig2) there. the distribution of the quantity below
    % should converge to that of the standard normal.
    % in the lecture slides, it should converge to N(0,sig2)
    sumX = (sum(X,2)-n*mu)/(sqrt(sig2*n));
    
    % we use 50 bins for our histogram
    numbins = 50;
   
    hist(sumX,numbins); %minval+stepsize/2:stepsize:maxval-stepsize/2);
    title (sprintf('PDF of the average of %d iid random variables',n));
    pause;
end
