%% central limit theorem
n = 10000;
m = 10;  % number of loops

%% (b) a single data set of random numbers:
x = rand(n, 1);

%% (c) plot probability density:
%histogram(x, 'Normalization', 'pdf');
[h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1));  % normalization
bar(b, h)
title('A uniform distribution')
xlabel('x')
ylabel('probability density')
pause(2.0)

%% (d) sum of two random numbers:
y = rand(n, 1);
x = x + y;
%histogram(x, 'Normalization', 'pdf');
[h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1));  % normalization
bar(b, h)
title('Sum of two uniform distributions')
xlabel('x')
ylabel('probability density')
pause(2.0)

%% (f) sum up more distributions:
x = zeros(n, 1);
means = zeros(m, 1);
stds = zeros(m, 1);
for i=1:m
    y = rand(n, 1);     % new uniform distributed numbers
    x = x + y;          % add them to the sum
    mu = mean(x);       % compute mean
    sd = std(x);        % compute standard deviation
    means(i) = mu;      % store mean
    stds(i) = sd;       % store standard deviation
    %xx = min(x):0.01:max(x);
    xx = -1:0.01:i+1;   % x-axis values for plot of pdf
    p = exp(-0.5*(xx-mu).^2/sd^2)/sqrt(2*pi*sd^2);  % pdf
    plot(xx, p, 'r', 'linewidth', 3 )
    ns = sprintf('N=%d', i);
    text(0.1, 0.9, ns, 'units', 'normalized')
    hold on
    %histogram(x, 20, 'Normalization', 'pdf');
    [h,b] = hist(x, 20);
    h = h/sum(h)/(b(2)-b(1));  % normalization
    bar(b, h)
    hold off
    xlim([-0.5, i+0.5])
    xlabel('x')
    ylabel('summed pdf')
    savefigpdf(gcf, sprintf('centrallimit-hist%02d.pdf', i), 6, 5);
    if i < 6
        pause(3.0)
    end
end

%% (h) mean and standard deviation in dependence on number of summed up distributions:
nx = 1:m;
plot(nx, means, 'b', 'linewidth', 4);
hold on
plot(nx, stds, 'r', 'linewidth', 4);
xx = 0:0.01:m;
sdu = 1.0/sqrt(12);   % standarad deviation of the uniform distribution
plot( xx, sqrt(xx)*sdu, 'k' )
legend('mean', 'std', 'theory')
xlabel('N')
hold off
savefigpdf(gcf, 'centrallimit-samples.pdf', 6, 5);