% Illustration of the central limit theorem (Chapter 4, Box B)
% Written by David Sumpter

%This code generates random numbers with various distributions and 
% compares their sum to the Normal distribution.

sizea=5000;
range=[-0.05:0.1:1.05];

x=(randn(sizea,25)+5)/10;
subplot(3,3,1)
[h tbin]=hist(x(:,1),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(x(:,1)),std(x(:,1)))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(a)')
subplot(3,3,2)
[h tbin]=hist(mean(x(:,1:5)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:5))),std(reshape(x(:,1:5),1,sizea*5))/sqrt(5))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(b)')
subplot(3,3,3)
[h tbin]=hist(mean(x(:,1:25)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:25))),std(reshape(x(:,1:25),1,sizea*25))/sqrt(25))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(c)')

x=binornd(1,1/2,sizea,25);
subplot(3,3,7)
[h tbin]=hist(x(:,1),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(x(:,1)),std(x(:,1)))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(g)')
subplot(3,3,8)
[h tbin]=hist(mean(x(:,1:5)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:5))),std(reshape(x(:,1:5),1,sizea*5))/sqrt(5))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(h)')
hlx=xlabel('Estimate')
set(hlx,'FontSize',14);
subplot(3,3,9)
[h tbin]=hist(mean(x(:,1:25)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:25))),std(reshape(x(:,1:25),1,sizea*25))/sqrt(25))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(i)')

x=(lognrnd(zeros(sizea,25)+1,ones(sizea,25)))/10;
subplot(3,3,4)
[h tbin]=hist(x(:,1),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(x(:,1)),std(x(:,1)))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
hly=ylabel('Frequency')
set(hly,'FontSize',14);
title('(d)')
subplot(3,3,5)
[h tbin]=hist(mean(x(:,1:5)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:5))),std(reshape(x(:,1:5),1,sizea*5))/sqrt(5))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(e)')
subplot(3,3,6)
[h tbin]=hist(mean(x(:,1:25)'),[0:0.1:1])
hold off
bar(tbin,h/sizea,'w')
nd=normcdf(range,mean(mean(x(:,1:25))),std(reshape(x(:,1:25),1,sizea*25))/sqrt(25))
nd2=nd(2:12)-nd(1:11);
hold on
plot(range(1:11)+0.05,nd2,'k')
axis([-0.1 1.1 0 1])
title('(f)')