%% Temple university CIS 2033
% A modern Introduction to Probability and Statistics
% Chapter 7: Expectation and variance
%
% Author: Djordje Gligorijevic
% email: gligorijevic@temple.edu
%
% Special thanks to Min Xiao (http://astro.temple.edu/~tud09663/index.html)
% for given examples and ideas for lab examples.

distnames = {'Uniform', 'GeneralizedPareto'};

%% For X ~ U(0, 1) compute true and empirical values of mean and variance

alpha1 = 0;
beta1 = 1;

% Compute TRUE values;
uds = makedist(distnames{1},'lower',alpha1,'upper',beta1);
mean1 = uds.mean;
variance1 = uds.var;

% Compute empirical values
repetition = 100; % repetitioneat 100 times;
sizes = [5,10,50,100,500,1000,5000]; % For different sampling sizes;
nS = length(sizes);
mean2 = zeros(nS,1);
variance2 = zeros(nS,1);
for sizeIndex = 1:nS
    aSize = sizes(sizeIndex);
    aMean = 0;
    aVariance = 0;
    for arepetition = 1:repetition
        ys = random(uds,[aSize,1]);
        aMean = aMean + mean(ys);
        aVariance = aVariance + var(ys);
    end
    aMean = aMean/repetition;
    aVariance = aVariance/repetition;
    mean2(sizeIndex) = aMean;
    variance2(sizeIndex) = aVariance;
end

fig_1a = figure;
xlim([min(sizes) max(sizes)]);
plot(sizes,mean2, '-r');
hold on;
plot(sizes,repmat(mean1,1,nS),'-b');
legend('EMPIRICAL','TRUE');
title('Mean');
saveas(fig_1a,'uni_mean.png','png');

fig_1b = figure;
xlim([min(sizes) max(sizes)]);
plot(sizes,variance2,'-r');
hold on;
plot(sizes,repmat(variance1,1,nS),'-b');
legend('EMPIRICAL','TRUE');
title('Variance');
saveas(fig_1b,'uni_var.png','png');

%% For X ~ Par(3) compute true and empirical values of mean and variance

alpha1 = 3;
k1 = 1/alpha1;
sigma1 = 1/alpha1;
theta1 = 1;

% Compute TRUE values;
uds = makedist(distnames{2},'k',k1, 'sigma',sigma1,'theta',theta1);
mean1 = uds.mean;
variance1 = uds.var;

% Compute empirical values;
repetition = 100;
sizes = [5,10,50,100,500,1000,5000];
nS = length(sizes);
mean2 = zeros(nS,1);
variance2 = zeros(nS,1);
for sizeIndex = 1:nS
    aSize = sizes(sizeIndex);
    aMean = 0;
    aVariance = 0;
    for arepetition = 1:repetition
        ys = random(uds,[aSize,1]);
        aMean = aMean + mean(ys);
        aVariance = aVariance + var(ys);
    end
    aMean = aMean/repetition;
    aVariance = aVariance/repetition;
    mean2(sizeIndex) = aMean;
    variance2(sizeIndex) = aVariance;
end

fig_2a = figure;
xlim([min(sizes) max(sizes)]);
plot(sizes,mean2, '-r');
hold on;
plot(sizes,repmat(mean1,1,nS),'-b');
legend('EMPIRICAL','TRUE');
title('Mean');
saveas(fig_2a,'par_mean.png','png');

fig_2b = figure;
xlim([min(sizes) max(sizes)]);
plot(sizes,variance2,'-r');
hold on;
plot(sizes,repmat(variance1,1,nS),'-b');
legend('EMPIRICAL','TRUE');
title('Variance');
saveas(fig_2b,'par_var.png','png');