%%%
% CLT central limit theorem 
% some demos
% Tomas Svoboda, 3/1999, CMP Prague
%%%

clear all

data.np	     = 50000;	% number of points in data set
data.densnp  = 2000;	% number of points for density function
demo.histnum = 30;		% number of bars in histogram

for i = 1:3,
  data.x{i}		= rand(data.np,1);		% data generation
  data.var(i)	= std(data.x{i},1)^2;	% variance of the data (maximum likelihood)
  data.mean(i)	= mean(data.x{i});		% mean of the data
  data.f{i}		= ones(data.densnp,1);	% density function of uniform distribution
end

%%%
% plots for x1
%%%
figure(1)
xaxis = [0:1/(data.densnp-1):1];
clf
hold on
n = hist(data.x{1},demo.histnum);
bar([1/2*1/demo.histnum:1/demo.histnum:(demo.histnum-1/2)/demo.histnum],n./mean(n));
plot(xaxis,data.f{1},'-g')
title('x_1')
axis([0 3 0 1.5])
text(2,0.9,'Density function \it f','Color','g')

%%%
% plots for x1+x2
%%%
f = conv(data.f{1},data.f{2});	% convolution of density function
f = f./(2*mean(f));				% density normalization
xaxis =2*[1/2*1/demo.histnum:1/demo.histnum:(demo.histnum-1/2)/demo.histnum];
% approximation by Central Limit Theorem
fclt = normpdf(xaxis,sum(data.mean(1:2)),sqrt(sum(data.var(1:2))));
x = data.x{1}+data.x{2};
figure(2)
clf
hold on
n = hist(x,demo.histnum);  
bar(xaxis,n./(2*mean(n)));
plot(xaxis,fclt,'-r','LineWidth',2);
plot([0:1/(data.densnp-1):2],f,'-g')
title('x_1+x_2')
axis([0 3 0 1.5])
text(2,1,'approx. by CLT','Color','r')
text(2,0.9,'Convolution of \it f','Color','g')


%%%
% plots for x1+x2+x3
%%%
f = conv(f,data.f{3});
f = f./(3*mean(f));
xaxis =3*[1/2*1/demo.histnum:1/demo.histnum:(demo.histnum-1/2)/demo.histnum];
fclt = normpdf(xaxis,sum(data.mean(1:3)),sqrt(sum(data.var(1:3))));
x = data.x{1}+data.x{2}+data.x{3};
figure(3)
clf
hold on
n = hist(x,demo.histnum);  
bar(xaxis,n./(3*mean(n)));
plot(xaxis,fclt,'-r','LineWidth',2);
plot([0:1/(data.densnp-1):3],f,'-g')
title('x_1+x_2+x_3');
axis([0 3 0 1.5])
text(2,1,'approx. by CLT','Color','r')
text(2,0.9,'Convolution of \it f','Color','g')
