%% 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 );