% copyright 2015 Keith Chugg, USC
%
% Compute the sample mean/variance as a function of n for 
% several realizations and plot 
%
% Also, for a given n, o=plot a histogram of the sample mean and variance
% 

clear all
clc

% The values below are for the plots of sample m, var vs. N
N=300;
n=1:N;
N_realizations = 100;

N_hist = 200;  % Plot histogram of sample m. sigma^2 @ this value of N
N_realization_hist = 100000;  % use many realizations for the historgams


%%%%%%%%%%%%%%%%%%%%
% Comment out the random number generator code below for all but the 
% one you want to see.
%%%%%%%%%%%%%%%%%%%%%

mode = 2;

if mode == 0
    % GAUSSIAN (0,1)
    m = 0;
    sigma = 1;
    x=random('norm',m,sigma,N,N_realizations);
    z=random('norm',m,sigma,N_hist,N_realization_hist);
end

if mode == 1
    % EXPONENTIAL 
    lambda_inverse = 4;
    x=random('exp',lambda_inverse,N,N_realizations);
    z=random('exp',lambda_inverse,N_hist,N_realization_hist);
end

if mode == 2
    % uniform
    a = -sqrt(3);
    b = sqrt(3);
    x=random('Uniform',a,b,N,N_realizations);
    z=random('Uniform',a,b,N_hist,N_realization_hist);
end

if mode == 3
    % Bernoulli
    p =0.3;
    x=random('Binomial',1,p,N,N_realizations);
    z=random('Binomial',1,p,N_hist,N_realization_hist);
end

%%%%%%%%%
% Below is to generate sample mean and variance vs n plots

sample_mean = zeros(N,N_realizations);
sample_variance = zeros(N,N_realizations);

for i = 1:N
    sample_mean(i,:) = mean(x(1:i,:));
    sample_variance(i,:)= var(x(1:i,:));
end

% Plot the sample mean/variance as function of n
figure
subplot(2,1,1)
plot(n,sample_mean)
title('Sample Mean')

subplot(2,1,2)
plot(n,sample_variance)
title('Sample Variance')

% Now plots histgrams of the sample mean and sample variance at one N

figure
subplot(2,1,1)
hist(mean(z),100)
title('Sample Mean using N Realizations - Histogram')

subplot(2,1,2)
hist(var(z),100)
title('Sample Variance using N Realizations - Histogram')


