%% (a) generate random numbers:
n = 100000;
x=randn(n, 1);

for nsamples=[3 5 10 50]
    nsamples
    %% compute mean, standard deviation and t:
    nmeans = 10000;
    means = zeros( nmeans, 1 );
    sdevs = zeros( nmeans, 1 );
    students = zeros( nmeans, 1 );
    for i=1:nmeans
        sample = x(randi(n, nsamples, 1));
        means(i) = mean(sample);
        sdevs(i) = std(sample);
        students(i) = mean(sample)/std(sample)*sqrt(nsamples);
    end
    
    % Gaussian pdfs:
    msdev = std(means);
    tsdev = 1.0;
    dxg=0.01;
    xmax = 10.0;
    xmin = -xmax;
    xg = [xmin:dxg:xmax];
    pm = exp(-0.5*(xg/msdev).^2)/sqrt(2.0*pi)/msdev;
    pt = exp(-0.5*(xg/tsdev).^2)/sqrt(2.0*pi)/tsdev;
    
    %% plots
    subplot(1, 2, 1)
    bins = xmin:0.2:xmax;
    [h,b] = hist(means, bins);
    h = h/sum(h)/(b(2)-b(1));
    bar(b, h, 'facecolor', 'b', 'edgecolor', 'b')
    hold on
    plot(xg, pm, 'r', 'linewidth', 2)
    title( sprintf('sample size = %d', nsamples) );
    xlim( [-3, 3] );
    xlabel('Mean');
    ylabel('pdf');
    hold off;
    
    subplot(1, 2, 2)
    bins = xmin:0.5:xmax;
    [h,b] = hist(students, bins);
    h = h/sum(h)/(b(2)-b(1));
    bar(b, h, 'facecolor', 'b', 'edgecolor', 'b')
    hold on
    plot(xg, pt, 'r', 'linewidth', 2)
    title( sprintf('sample size = %d', nsamples) );
    xlim( [-8, 8] );
    xlabel('Student-t');
    ylabel('pdf');
    hold off;
    
    savefigpdf( gcf, sprintf('tdistribution-n%02d.pdf', nsamples), 14, 5 );
    pause( 3.0 )
end